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The  objective  of  the  study  reported  herein  was  to  develop  a  general  soil  stress-strain 
w" ’ r h  can  be  used  to  solve  a  wide  range  of  soil  dynamics  problems  of  interest  to  the 
'■ir  Force.  The  approach  used  was  to  review  existing  soil  constitutive  models  usecj  to  predict 
•re  response  of  soil  masses  to  complex  dynamic  loads,  and  then  formulate  a  new  model  for  that 
.jrpose.  Eight  existing  soil  dynamic  stress-strain  models  were  studied,  including  exercising 
•  >'<"  .L'lni  i  i "’ii". i n  stress  and  strain  paths  for  comparison.  The  models  were:  linear  elastic, 
i  f‘<i  r  v  i  si  of  hist  ii. ,  hyperbolic,  F’yke  cyclic  simple  shear,  elastic-perfectly  plastic,  modified 
•rWL  engineering,  effective  stress  cap,  and  Lade.  Each  model  is  discussed  and  reviewed.  The 
iiscussion  of  each  model  includes:  motivation,  assumptions,  basic  equations,  parameter 
letermination,  and  computed  behavior.  Based  on  the  above  review,  the  Lade  model  was  selected 
: s  the  best  point  of  departure  for  developing  a  new  soil  stress-strain  model  for  complex 
iynamic  loading,  because  of  its  accuracy  and  flexibility  in  representing  soil  stress-strain 
uehavior,  ease  of  parameter  determination,  and  ease  of  developing  intuition  for  parameter 
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ABSTRACT  (continued) 

physical  significance  and  accuracy.  The  same  five  attributes  discussed  for  each  existing 
model  were  examined  for  the  new  ARA  conic  model,  so  called  because  its  principal  mathemat¬ 
ical  surfaces  are  conic  sections.  The  computer  code  used  to  exercise  all  nine  soil  con¬ 
stitutive  models  under  eleven  stress  and  strain  paths  is  called  the  Soil  Element  Model 
(SEM).  It  can  be  incorporated  in  large  finite  difference  or  finite  element  codes  for 
analyzing  the  response  of  soil  masses  to  complex  dynamic  loads. 

The  ARA  conic  model  performs  well  over  a  wide  range  of  loading  conditions,  many  depar¬ 
ting  considerably  from  those  used  to  determine  the  model  parameters.  The  parameters  are 
determined  in  a  straightforward  manner,  and  the  model  reflects  the  influence  of  the  inter¬ 
mediate  principal  stress  on  shear  strength  through  a  shear  failure  surface  involving  three 
I  independent  stress  invariants:  the  first  total  stress  invariant  and  the  second  and  third 
Ideviator  stress  invariants.  For  this  reason  the  model  is  also  called  a  three  invariant 
model.  Measured  shear  strengths  in  both  compression  and  extension  can  be  matched  exactly, 
and  the  mathematical  formulation  of  the  shear  failure  surface  is  such  that  the  shear 
strength  for  any  value  of  the  intermediate  principal  stress  can  be  computed  directly 
[without  trial  and  error.  The  ARA  conic  model  also  exhibits  dilatancy,  generates  only 
'positive  plastic  work,  and  has  a  provision  for  strain  softening  in  shear. 
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(A. 4b) 


(A. 4ci 


Because  the  area  of  a  pyramid  cross  section  taken  parallel  to  the 
base  is  proportional  to  the  square  of  the  perpendicul ar  distance  to  the 
cross  section  from  the  apex,  the  volume  of  a  pyramid  is  given  by  the  well 
known  formul a: 

volume  =  ^e1 1?-1'*'  x  base  area 


Therefore,  the  volume  of  tetrahedron  OABC  can  be  expressed  in  any  of  four 
ways: 


V  =  (area  ABC)  =  (area  BOC)  =  (area  AOC) 
from  which  it  follows  that,  using  Equations  (A. 4), 


[area  AOB) 


area  BOC  h^  _ 

area  ABC  a  “in 


area  AOC 


area  AOB  h  / «  t 

area  ABC  c  a3n 

Consider  now  the  stresses  exerted  on  tetrahedron  OABC  by  the 
surrounding  material,  shown  in  Figure  (A. 2).  The  stress  vectors  on> 

-o^,  -o^  and  -o^  are  equal  to  the  resultant  forces  acting  on  faces  ABC, 
BOC,  AOC  and  AOB,  respectively,  divided  by  the  area  of  the  face  on  which 
they  act.  In  the  limit,  as  the  dimensions  of  the  tetrahedron  approach 
zero,  the  stress  vectors  each  represent  a  uniformly  distributed  force  per 


unit  area.  Equilibrium  considerations  dictate  that  the  vector  sum  of  all 
forces  acting  on  the  tetrahedron  be  zero.  Body  forces,  such  as  those  due 
to  gravity  and  inertia,  need  not  be  included  in  the  vector  equilibrium 
equation  for  an  infinitesimal  tetrahedron,  since  their  magnitudes  are 
proportional  to  h  ,  whereas  the  magnitudes  of  surface  forces  are 

p 

proportional  to  h  .  The  body  forces  therefore  approach  zero  more 
rapidly  than  do  the  surface  forces  as  the  tetrahedron  dimensions  approach 
zero.  Therefore,  the  following  development  is  as  applicable  to  problems 
of  wave  propagation  as  it  is  to  static  problems.  The  vector  equilibrium 
equation  for  the  infinitesimal  tetrahedron  is: 

~n  (area  ABC)  -  (area  BOC)  -  03  (area  AOC)  -  03  (area  AOB)  =  0  (A. 7 

Dividing  both  sides  of  Equation  (A. 7)  by  (area  ABC)  and  using 
Equations  (A. 6)  yields 

on  =  Oid.  +  oQa9„  +  a-,ao„  =  o.a. 

n  I  In  2  2n  j  Jn  j  jn 

Equation  (A. 8)  is  commonly  referred  to  as  Cauchy's  equation. 

Each  of  the  stress  vectors  in  Figure  (A. 2)  can  be  written  in 


component  form,  so  that 


°ij  =  stress  acting  in  direction  ei ,  on  the  face  with  outward 


uni t  normal  e • . 

J 

Equations  (A. 9)  can  be  written  concisely  in  the  form: 

n  .  =  a  ■  .e  . 

J  1  J  1 

Substitution  of  Equation  (A. 10)  into  Equation  (A. 8)  yields: 


(A.lt 


=  o -a  .  =  (a  .  .e  .  )a  • 

J  jn  i j  i  jn 


=  ( o  •  •  a  •  )  e  • 
ij  jn'  l 


or,  in  matrix  form: 


( A.  1 


{°n)  =  £  {anl 

The  component  of  on  in  the  direction  of  the  unit  vector  m  in 
Figure  (A. 2)  is: 

°mn  =  m  an  =  {ap}T  {on}  =  {am)T  a{an)  =  aimoijajn 


(A.i; 


( A.  1 


Consider  now  two  rectangular  Cartesian  coordinate  systems  having  a 
common  origin,  as  shown  in  Figure  (A. 3),  and  let 

a..  =  stress  component  defined  in  the  1,2,3  coordinate  system 

*  J 

o.!j  =  stress  component  defined  in  the  r,2',3’  coordinate  system 
a.  .  =  cosine  of  the  angle  between  e  .and  e' . 

*  *j  '  0 

Then  Equation  (A. 13)  gives 


,  T 

u  =  a  .  o  .  .  a  .  -  a  .  o  •  •  a  • 

Rin  im  ij  jn  rrn  ij  jn 
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From  Equation  (A. 13)  it  also  follows  that  the  normal  component  of  an  in 
Figure  (A. 2),  denoted  inn>  is  simply 
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nn 


=  n 
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a  r  a  i 

— i  n  • 


(A. 16) 


Principal  Stresses 

It  is  of  interest  to  investigate  the  change  in  .Tnn  in 
Equation  (A. 16)  due  to  an  infinitesimal  change  i n  77 .  Figure  (A. 4)  shows 
that  since  77  must  remain  a  unit  vector,  it  can  only  change  direction,  so 
that 

/da  ;  =  dXj-a^i  (A.  17) 

n  J  it. 

where 


rv 


:an) 


=  0 


(A. 18) 


Equations  (A. 17)  and  A. 18)  merely  state  that  an  infinitesimal  change  in 
the  unit  vector  77  must  lie  in  the  plane  normal  to  77.  The  increment  in 
7nn  due  to  dn  can  now  be  obtained  by  differentiating  Equation  (A. 16), 
which  yi el ds 


da 


nn 


=  {darTr  -l!an )  +  ran|T  -2.{dan7  =  ldan}T  +  {anl  (A. 19) 


However,  because  there  are  assumed  to  be  no  distributed  torques,  the 
stress  matrix  is  symmetric,  i.e.: 

aT  =  a  (A. 20) 


Therefore,  substitution  of  Equations  (A. 17)  and  (A. 20)  into 
Equation  (A. 19),  taking  account  of  Equation  (A. 13),  yields 


in  the  direction  of  t.  From 


where  otn  is  the  shearing  component  of  on 
Equation  (A. 21)  it  follows  that  stationary  values  of  normal  stress  occur 
on  plane?  of  zero  shear.  Therefore,  to  locate  stationary  values  of  normal 
stress,  one  seeks  those  directions  for  which  the  stress  vector,  an,  has 
no  component  normal  to  n,  and  is  therefore  parallel  to  n.  For  such  a 
"principal"  direction, 


°n  =  cn 


(A. 22) 


Using  Equation  (A. 12),  Equation  (A. 22)  can  be  written  in  matrix  form  as 


{V 


=  £{anj  -  °{V 


(A. 23) 


or,  in  homogeneous  form 

(c  -  an|  =  { o } 


(A. 24) 


The  expanded  form  of  Equation  (A. 24)  is 


0 11  -  0 

°12 

°13 

f°ln] 

°2 1 
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V2nf 
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°33  '  a 

(a3n) 

(A. 25) 


and  since  Equations  (A. 25)  are  a  set  of  homogeneous  linear  equations,  they 
possess  a  unique  solution  if  and  only  if  the  determinant  of  the 
coefficient  matrix,  in  this  case  referred  to  as  the  character!' sti c 
determinant,  is  zero.  The  equation  expressing  this  condition  is  the 
characteristic  equation: 


Expansion  of  Equation  (A. 26)  yields 


or 

where 
equati on. 
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3  2 

-(n  -  I  jo  -  l^rj  -  1^)  =  -(o  -  o^)(a  -  o^Ho  -  a^)  =  0  (A. 2 

o ^ ,  and  are  the  three  real  roots  of  the  cubic  characteristic 

called  the  principal  stresses,  and  the  stress  invariants,  I j, , 


lo  and  U  are  defined  as  follows: 


°21  °22 

°23 

°31  °32 

c33 
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Cardan's  classic  solution  of  Equation  (A. 27)  [USPENSKY  (1948:84)]  begins 
by  setti ng 


o  =  s  + 
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(A. 


which  causes  Equation  (A. 26)  to  take  the  form 


or 
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=  0 


(A. 


where  the  deviator  stress  matrix,  s,  is  defined  by  the  equation 


Expansion  of  Equation  (A. 33)  yields 


(s  -  J ^ s  -  J3)  =  -  (s  -  s^Hs  -  s2)(s  -  S3)  =  0 


(A. 3 


where  s^,  s2,  S3  are  the  three  real  roots  (principal  deviator  stresses), 
and  the  deviator  stress  invariants,  ,  J2,  and  J3  are  defined  as  follows 

(A.3i 


+  s22  +  s33  = 

S1  +  s2  + 
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The  s2  term  is  absent  from  Equation  (A. 35)  because  is  zero,  and 
Equation  (A. 35)  is  therefore  referred  to  as  the  reduced  characteristic  or 
reduced  cubic  equation. 

The  solution  to  Equation  (A. 35)  is  assumed  to  be  of  the  form 

s  =  A  +  B  (A. 3! 

so  that 


s3  =  A3  +  3A28  +  3AB2  +  B3 


A3  +  B3  +  3AB( A  +  B) 


A3  ♦  B3  +  3ABs 


or 


s3  -  3ABs  -  (A3  +  B3)  =  0 


( A .  4( 


Comparison  of  Equations  (A. 40)  and  (A. 35)  shows  that 


3AB  =  J. 


U 


or 


A3B3  = 


T~ 


U 


and 


3  3 

AJ  +  B-3  =  J- 


(; 


Thus,  the  sum  and  product  of  A3  and  B3  are  specified,  which  means  that 
A3  and  B3  are  the  roots  of  a  quadratic  equation.  To  obtain  the 
quadratic  equation  we  first  obtain  an  expression  for  B3  from 
Equation  (A. 43) . 

(/ 


3  3 

B3  =  J3  -  R 


Substitution  of  Equation  (A. 44)  into  Equation  (A. 42)  then  yields 


A3(J  -  A3)  =  J/3  -  (A3)2  =  -J- 
3  3  3 


or 


,  .3.2  .  .3  .  J2| 

(A  }  -  J3A  \T 


Because  of  the  symmetry  of  Equations  (A. 42)  and  (A. 43),  B3  also 
satisfies  Equation  (A. 45).  Therefore,  the  solutions  for  A3  and  B3  are 


3  °3 
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~  TT 


\T 


106 
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In  addition,  Equation  (A. 36)  requires  that 

Si  +  s2  +  s3  =  3r  ♦  uq  =  o  (A.! 

so  that  the  deviator  stress  paramenters  r  and  q  are  related  by  the 
expression 

r  =  •  j  q  (A.! 


and  therefore  the  principal  deviator  stresses  can  all  be  expressed  in 
terms  of  the  parameters  u  and  q  as  follows: 


1  =  (1  - 

$><1 

(A. 

2  =  7^ 

(A. 

3  "  -  (1 

+  ^ 

(A. 

Substitution  of  Equation  (A. 56-58)  into  the  ratio  in  question  reduces  it 
to  a  function  of  u  only,  as  follows: 


According  to  Equation  (A. 50)  the  parameter  4  varies  between  -1  (when 

s2  =  S3)  and  +1  (when  s2  =  s^).  The  ratio  in  Equation  (A. 61) 
therefore  also  varies  between  +1  and  -1,  as  can  be  verified  by  direct 
calculation,  and  therefore  the  radicand  in  Equations  (A. 46)  is 
nonpositive.  This  being  the  case,  the  quantities  A^  and  are 
complex,  and  can  therefore  be  represented  by  an  Argand  diagram  as  shown  in 
Fi gure  ( A. 5) . 

The  expressions  for  AJ  and  B  ,  shown  in  Figure  (A. 5)  are: 
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(A. 62a) 
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(A. 62b) 


where 


cos  3u  = 


2) 

\Ti 
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:  A.63) 


(3 
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Equation  (A.63)  is  plotted  in  Figure  A. 6. 

The  desired  values  of  A  and  B  are  therefore  the  complex  cube  roots  of 


the  expressions  in  Equations  (A. 62),  i.e.: 


Now  Equation  (A. 39)  states  that  each  of  the  three  roots,  s^  s^,  and 
s^  is  the  sum  of  two  quantities,  A  and  D,  whose  product  is  real  according 
to  Equation  (A. 41).  Therefore,  the  values  of  A  and  B  used  to  form 
s.  (j  =  1,2,3)  must  be  complex  conjugates,  which  requires  that  k  =  j  in 

J 

Equations  (A. 64).  Equation  (A. 39)  therefore  takes  the  form: 

J  P  i  u  j  -  i  w  •  J  p 

s.  =  l-j-  (e  J  +  e  J)  =  cos  <jj  (j  =  1,2,3)  (A. 6 


or 
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( A .  61 


(A. 6 1 
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s,  =  2 


COS 


so  that  Equation  (A. 31)  takes  the  form: 


*1  ^2  *1 

S1  +  T~  =  2Jl~  cos  ,J1  +  T 


h  J2  h 

°2  =  S2  +  T  =  2Jt~  C0S  iJ2  +  T 


^1  r2  *1 

°3  =  s3  +  y  =  2Jl~  cos  'J3  +  T 


(A. 69 


(A. 69 


(A. 69 


principal  Stress  Directions 

Once  the  principal  stresses  have  been  calculated,  the  direction 
cosines  associated  with  each  principal  stress  can  be  obtained  from  three 
separate  solutions  of  Equation  (A. 25)  in  which  a  is  successively  set  equal 
to  uj,  o2  and  o ^ ■  In  each  solution  at  least  one  of  the  equations  will  be 
redundant,  so  that  a  corresponding  number  of  the  direction  cosines  will  be 
arbitrary.  The  directions  thus  obtained  are  those  of  the  principal  stress 


axes,  and  can  be  shown  to  be  orthogonal  as  follows.  Let 


2-  =  \{al}  \  {“2}  |  {a3)J 

where  ja^j  ,  |a2|  and  jc^j  are  the  three  principal  directions.  Then 
Equation  (A. 23)  can  be  written  in  the  form: 


(A. 70 


Premultiplying  both  sides  of  Equation  (A. 71)  by  a  ,  then 
transposing  the  result  and  taking  account  of  Equation  (A. 20)  yields 
T  T  r 

a  o  a  =  a  a  ■  o 

-  - - Pj 


T  T  T 


T  r  T 

a  o  a  =  fCp^  a  a 


a)  ("a  =  G  (a"'"  a) 

_  _  ^  Pj  -  - 


The  only  way  Equation  (A. 73)  can  be  satisfied  is  for  the  matrix 

a  to  be  diagonal,  which  means  the  columns  of  £  are  orthogonal 

(perpendicular).  Since  the  {aj }  are  unit  vectors,  we  have 

T  T 

a  a  =  1 


so  that  Equation  (A. 71)  yields 


APPENDIX  B 


CAYLEY-HAMILTON  INVARIANT  FORMULATIONS 


If  a  is  the  3  x  3  matrix  containing  the  column  vectors  of  direction 

cosines  of  unit  vectors  in  the  three  principal  stress  directions  of  the 

stress  matrix,  o,  and  rr,  1S  the  diagonal  matrix  of  corresponding 

u 

principal  stresses,  then  Equations  (A. 74-76)  give 


T 

a 


a 


=  I 


T 

a  a  a  = 


ra 


Pj 


It  is  easy  to  generalize  Equation  (A. 76)  for  higher  powers  of 
ana  a.  For  example, 


(A. 

(A. 

(A. 

CP 


r-  2  T  ,  _  I,,  r  T,  _  2 

a  *0,.  a  =  (a  a  )(a  la  a  )  -  a_ 

-Pj-  -  Pj-“Pj 


and ,  i n  general , 


a  Hs  n  aT  =  (a  HI  aT )  =  a”  (n  an  integer)  (B. 

-  Pj  -  -  Pj- 

Now  if  Equation  (A. 27)  is  written  for  each  of  the  three  principal 
stresses,  and  the  result  put  in  diagonal  matrix  form,  we  obtain 

r„  3  -  I.  rn  2  -  I  n  -  Iol  *  ro,  (B. 

Pj  1  ^  2  pj  3-  -t 

Premultiplication  of  Equation  (B.3)  by  a  and  postmultiplication  by 


then  yiel ds 


(B.4) 


£3  -  Ii£2  -  ^£  -  1^2  =  2 


Equation  (B.4)  states  that  the  stress  matrix,  £,  satisfies  the  matrix 
form  of  its  own  characteristic  equation.  This  is  a  particular  example  of 
the  Cayley-Hamil ton  theorem,  and  leads  to  useful  expressions  for  the 
invariants  Ij,  l^,  I^,  Jj,  and  J^,  as  shown  below. 

Recall  that  the  inverse  of  a  square  matrix  equals  its  adjoint  (the 
transpose  of  its  matrix  of  cofactors),  divided  by  its  determinant.  Thus, 

vT 

—  (B.5) 


■1  r 


Now 


1 2  =  -  Tr(£)  =  -  Tr(T') 


(B.6 ) 


and 


:3  =  £l 


(B.7) 


so  that 


Tr(o_1)  = 


T r(  IT) 


2 

r3 


(B.8) 


Now  first  note  that 


Ij  =  Tr(o) 


(B.9) 


Next,  premultiply  Equation  (B.4)  by  £-1  to  obtain 

-  ~  !1-  '  Z2-  "  I^-"1  =  -  (B.10) 

Taking  the  trace  of  Equation  (B.10)  then  yields 
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Finally,  taking  the  trace  of  Equation  { B . 4 )  yields 
Tr(^3)  -  I1Tr(a2)  -  -  3I3 

=  Tr(a3)  -  IJrU2)  -  -^LTr(a2)  -  I2]  -  3I3 

l3 

=  Tr(a3)  -  lijrfa2)  +  2^  -  3I3  =  0  (B 

so  that 

I3  =  i[2Tr(o3)  -  SIjTrU2)  +  I3]  (B 

Since  the  deviator  stress  matrix,  s^  can  be  viewed  as  a  particular 
type  of  stress  matrix  whose  trace  {first  invariant)  is  zero. 

Equations  (B.9),  (B.12),  and  (B. 14)  yield 


For  machine  computation,  the  following  invariant  expressions  are 
convenient  because  they  do  not  require  calculation  of  principal  stresses: 


ll  =  °11  +  °22  +  °33 


(s  2  +  S  2  +  s  2)  +  2(s  2  +  s  2  +  s  2) 
($H  s22  s33  ;  s23  $3i 


(B. 16) 


(B.19) 


(8.20) 


3  =  11  22  33  +  ^12  23  31  +  s32s21b13 

2  2  2 
-  <  sns23  +  s22s31  +  S 33 S 1 2  1 

Equation  ( B .18)  is  identical  to  Equation  (A. 28);  Equation  (B.19)  is 
an  expansion  of  Equation  (B.16);  and  Equation  (B.20)  is  an  expansion  of 
Equation  (A.38). 

Another  useful  expression  for  J2  is  obtained  by  setting 


J2  -  J2  - 


2  2  2  2 
3(sll  +  s22  +  s33}  ‘  (sll  +  S22  +  S33}  .  ,2  +  c2  +  c2  * 

- 5 -  s12  s23  s31; 


2  2  2 

2 < s 1 1  +  s22  +  s33  ‘  S11S22  "  S22S33  '  S33sll)  +  ,2  +  2  +  2  , 
- 5 -  s12  s23  S31 ' 


2  2  2 

(sll  -  s22)  +  (s22  '  s33^  +  (s33  "  Sll)  +  /r2  +  c2  t  t2  , 
- 5 -  s12  s23  S3r 

(o  -  o22)2  +  (o22  -  o33)  +  (o33  -  on)  2  2  2 

- Z -  °12  23  °31J  lB'21) 


By  using  Equation  (A. 34)  in  the  form 


Equations  { B . 1 2 )  and  (B.14)  can  be  made  to  yield 


6' 


and 


1 2  =  yLTr(o2)  -  I2j  =  ^Jlr[(s  ♦  ^  U2]  -  I2J 


lr 

7 


2  ^  21 1 


[T r(_st  +  —y-  s  +  -g— _!_)  -  I ,  j 


r2, 

'  1  ' 


=  J2-~ 


(B.2 


1 3  =  y[2Tr(o3)  -  3 1 L T r ( o_2 )  ♦  I3] 


■g,|2Tr[(_s  +  y-  _I_)  j  -  3IjTr[(  +  y-  _I_)  J  +  Ij 


I2  I3 

A  1  1  - 


=  i[2Tr(s3  +  Ixs2  +  yU  ♦  ^rl)  -  3I1Tr(s2  +  +  +  ij] 


2I1  ll  J1  *1 

=  J3  +  "rt  +  77  ’  IiJ2  '  r-  +  r* 


!1J2  !1 
J3  ■  ~3  +  77 


(B.2' 


so  that 


APPENDIX  C 


OCTAHEDRAL  PLANE  PLOTS 

The  positive  directions  of  the  orthogonal  principal  stress  axes, 
defined  hy  Equation  (A.70),  are  arbitrary.  Thus,  the  principal  stress 
axes  define  eight  octants  within  which  the  stress  vector  repeats  itself  in 
a  symmetrical  fashion.  It  is  convenient  to  use  the  principal  stress  axes 
as  coordinate  axes,  as  shown  in  Figure  (C.l).  The  principal  stress  axes 
define  the  principal  stress  space,  and  in  that  coordinate  system  the 
strpss  matrix  reduces  to  the  diagonal  matrix  of  principal  stresses. 


al 

0 

0 


(C.l) 


The  hydrostatic  axis  in  principal  stress  space  has  the  equation 


al  -  n2  ~  a3 


(C.2) 


and  the  unit  vector  pointing  away  from  the  origin  along  the  hydrostatic 
axis,  having  equal  positive  components  along  each  of  the  three  principal 
stress  axes,  is  the  octahedral  unit  normal. 


n0CT 


+  e35 


(C.3) 


The  plane  through  a  stressed  point  for  which  nggy  is  the  unit 
normal  vector  is  called  the  oct.rahpdral  plane.  Referring  to 
Equations  (A. 8)  or  (A. 12),  the  octahedral  stress  vector,  i.e.,  the  stress 
vector  acting  on  the  octahedral  plane,  is 


The  noma!  component  of  the  octahedral  stress  vector,  called  the 


octahedral  normal  stress,  is 


a0CT  "  a0CT  •  n0CT 


"l  +  ♦  a3  !1 


and  the  tangential  or  in-plane  component  of  oggy,  called  the  octahedral 
shear  stress  vector,  is 


t0CT  ■  °0CT  "  (Vt'  n0CTl  n0CT 


i  —  —  —  °nr  t  —  —  — 

=  yFr  (alel  +  n?e2  +  a3e3)  "  yy  (el  +  e?  +  e3 ^ 


yf  (Slel  '  4  S3e3> 

The  magnitude  of  the  octahedral  shear  stress  vector,  called  the 


octahedral  shear  stress,  is 


t0CT  =  J  t0CT  t0CT 


7  2  7  - 

s,  +  si  +  si  2J, 

1  t  3  c 


so  that  Equations  (A. 68)  can  be  written  in  the  form 


•i  - 


Tqct  cos  «j 


toct  cos  u2 


Tqct  cos  «3 


When  plotting  triaxial  stress  data  it  is  convenient  to  work  with  a 
mathematical  vector,  ”,  called  the  principal  stress  vector,  shown  in 
Figure  (C.2),  the  components  of  which  are  the  principal  stresses,  i.e. 


:7  -  rr  p  +  />r  p  +  a  p 

OCT  '  11  22  33 


(C.9) 


It  is  also  convenient  to  define  a  second  set  of  mutually  orthogonal, 
right-handed  unit  vectors,  called  the  octahedral  unit  vectors,  according 
to  the  relations 


n3  '  n0CT 


'  el  *  e2  *  e3) 


nl  *  jy  1  el  "  e3* 


(C.10) 


(C.ll) 


Vn3xnr 


1  0  -1 


el  *  ?e2  -  e3 


(C.12) 


The  unit  vectors  n^  and  n p  are  both  normal  to  n^y  ^nd  therefore 
both  lie  in  the  octahedral  plane.  The  unit  vector  has  a  positive  ej 
component,  an  equal  and  opposite  e^  component,  and  no  e?  component.  The 
nj,  rip,  and  n^  components  of  the  principal  stress  vector  are 


i  i 


X  =  r,  .  n1  = 


nl  "  "3  sl  *  s3 


(C. 13) 


y  =  n  •  P.p  = 


2°2  "  °1  "  a3  2s2  -  S1  -  5 


1  '  s3  fT 

=  h  ' 


(C.  14) 


z  =  a  .  n-j  = 


al  +  a2  +  a3 


» 


(C. 15) 


■  -  d 

!  -A 


r»r 


>.  WK  'U'MIifM  ■■  B'»  .■  .■  Km  J  V 


I 


I 


The  octahedral  plane  component  of  the  principal  stress  vector  is 


R=CT-  (n  •  n^) 


1^1  ,J?C?  °3e?)  '  ft  r0CT  ‘  JY  ^ 


(",e,  +  a„e„  + 


,e]  +  e?  *  e3) 


"  Slel  +  S2e2  +  s3e3 


and  the  magnitude  of  the  octahedral  plane  component  of  the  principal 
stress  vector  is 


(C . 16 ) 


R  =  jrrj  .  M  +  s?2  +  s?3  =fir2  =  jj 


OCT 


(C.17) 


It  can  easily  he  verified  by  direct  expansion,  using  Equations  (C.17) 
and  (C.l*),  that 

2  >  ? 

3  p  s,  -  2SjS3  +  3s? 

x  +  y  _  +  — 


s^  -  2s  s  +  s^  +  2s^  +  (-s  -s 
S1  ^sls3  3  2  1  S1  s3 


s"  +  s"  +  s"  =  RL  =  2J0  =  3  t2 


;2  ♦  S2  ♦  S2 
’1  s2  s3 


OCT 


(C. 18) 


Nov/  Equations  (C.13)  and  (C.14)  also  yield 


y  l|!"rV  1  l?h  '  %1  '  =3 


■/TV  ”1  -  "3  j  /T  V  si  -  s3  J  JT 

where  u  is  Lode's  parameter,  defined  by  Equation  (A. 50),  so  that, 
referring  to  Figure  (C.3), 

y  =  R  (—* — \ 
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(C. 1°) 


(C.20) 


*  ^V*V*jvV,V*  .  •* 


• 


*  -■  ' 


v;  ■"« 


■  .•  1 

■  i 


u  ■*-  ■'-  -'•-  ■'-  ■'  V 


or 


t 


i 


i 


or 


OCT 


1/3  +  y 


Comparison  of  Equations  (C.8b)  and  (C.21)  shows  that 


COS  Up  “ 


3  +  y‘ 


which  can  be  checked,  since 


cos  3u£  =  4  cos  u2  -  3  cos  <j2 


3 


4y"  3U 

r??77 '  TT^Tj177 


(3  +  y") 


4U3  -  3y  (3  +  y2)  y3  -  % 

- 7 7377 


(3  +  y2)3/2 


(3  ♦  4  r 


(3  +  y ) (3  -  y) 


(C.21 


(C. 


(C.23 


which  agrees  with  Equation  (A. 63).  Thus,  Equations  (C.13)  and  (C.14)  can 
be  written  in  the  form 

x  =  ft  rQCT  sin  u2  (C.24 


=  J7 


t0CT  C0S  'j2 
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y 


(C.2 


f\> 


A  common  equation  for  soil  shear  strength  is  the  Mohr- Coulomb 
equation,  which  can  be  written  in  the  form 

q=(d  +  p)sini6  ( 
where  the  friction  angle,  0,  may  depend  on  u;  d  is  a  constant  material 
parameter;  and 


Now  Equations  (C.13)  and  (C.14)  can  also  be  written  in  the  form 

x  =  /(T  q  ( 

2(o.  a0  .  o-,l  -  3(o,  -  o^)  . — 

>  *  — - 1 - J - 1 - L  =  f  <°0CT  -  Pi  I 

Substituting  Equation  (C.26)  into  Equation  (C.29),  then  solving  for  p 


y i el  ds 


x  =  (d  +  p)  si n  (5 

p  =  r  *---  d 
-V  2  sin  i 


ana  substituting  Equation  (C.32)  into  Equation  (C.30)  yields 


y  =  r  °oci 


—  'fa  (  d  *  O  Qp  T  )  - 


sin  i> 


Fi nal ly ,  if  we  set 


then  Equations  ( C . 1 9 )  ana  tC.33)  yield 


li  =  u  (C.3 

x 1  JT 

y  '  =  {t~ - —  x'  (C.3 

sin? 

hr'.trt,  when  a  =  U, 


air  i 


(C.3 

.Vi  ,how  a  rapid  method  for  plotting  triaxial 
;r  3  h„'t,  OaSea  on  Equations  (C.36)  and  (C.37) 

equation  it. 36)  applies  [Merkle  (1971:346)] 
■:  •_  ,jl  jficar  failure  criterion  than  the 
•<  ir  which  the  octahedral  shear  stress  is  a 
'  orr.iu  1  stress  and  the  octahedral  polar  angle, 


’  •  •  '  (C.3 

)  r,  this  Last-  a  unique  oetaheoral  cross  section  can  be  obtained  by 


I  fV 


APPENDIX  D 


BASIC  EOUATIONS  OF  ELASTOPLASTICITY 


If  a  cylindrical  soil  specimen  is  consolidated  under  an  isotropic 
stress  (ajc  =  =  a  ),  then  subjected  to  drained  compressive 

loading,  unloading,  and  reloading  under  constant  confining  stress 

=  ~3  =  ~3c’  ~1  stress-strain  curve  appears  as 

shown  in  Figure  (D.l).  Several  important  features  are  shown  in 
Figure  (D.l): 

1.  The  stress-strain  curve  is  nonlinear,  even  for  small  stresses  and 
strai ns. 

2.  Upon  unloading  from  point  A,  some  of  the  total  strain  is 
recoverable  (BE),  but  the  remainder  is  irrecoverable  (OB). 

3.  Reloading  occurs  along  a  path  (BC)  somewhat  different  from  the 
unloading  path  (AB),  until  nearing  the  previous  maximum  stress. 

At  this  point  additional  loading  approaches  and  proceeds  along 
what  appears  to  be  a  continuation  of  the  virgin  compression  curve 
(OA),  with  little  apparent  further  influence  of  previous 
unloading  or  reloading. 

4.  Unloading  and  reloading  occur  along  paths  whose  secant  from  zero 
to  maximum  stress  has  a  slope  very  close  to  that  of  the  initial 
tangent  to  the  stress-strai n  curve.  This  means  that  the 

i recoverable  portion  of  any  strain  increment  is  essentially  the 
difference  between  the  total  strain  increment  and  the  strain 
increment  associated  with  a  straight  line  stress-strain  curve 
having  a  slope  equal  to  that  of  the  initial  tangent  to  the  actual 
stress-strain  curve. 

By  convention,  recoverable  strains  are  called  elastic,  and 
i recoverable  strains  are  called  plastic.  If  the  unloading  and  reloading 
curves  in  Figure  (D.l)  both  retraced  the  virgin  loading  curve  (OA)  instead 
of  following  curves  (AB)  and  (BC)  the  stress-strai n  behavior  would  be 
called  non! i nearly  elastic.  As  it  is,  the  linear  portion  of  the  stress- 
strain  behavior  shown  in  Figure  (D.l)  is  elastic,  and  the  nonlinear 
portion  is  P1astic.  Since  some  of  the  stress-strai n  behavior  is  elastic 


and  the  rest  is  plastic,  the  overall  stress-strain  behavior  is  called 
elastoplastic. 


Multiaxial  elastoplasticity  theory  extrapolates  the  above  one- 
dimensional  stress-strain  observations,  arid  assumes  that  plastic  strains 
are  superimposed  on  elastic  strains.  Thus 

e.  .  =  £e-  +  CP.  (D.l) 

lj  lj  U 

with  a  similar  relation  holding  for  each  strain  increment. 


de . .  =  dcT,  +  dt; ; 


1J 


1  j 


1  J 


(D.2) 


When  the  elastic  behavior  is  isotropic,  the  stress  increments  are 
related  to  the  elastic  strain  increments  by  the  equations 

da. .  =  C®  ,  dcj,  (D.3) 

ij  ljkl  kl 


where  the  elastic  incremental  stiffness  tensor  is  given  by  the  expression 


cTjki  *  "ViAl  *  M»-Ko> 


(D.4) 


and  M  =  constrained  elastic  modulus 

K0  =  coefficient  of  elastic  lateral  stress  at  rest. 

The  parameters  M  and  KQ  are  sometimes  assumed  to  be  constant,  independent 
of  strain,  and  are  sometimes  assumed  to  vary  in  a  prescribed  manner. 

Once  the  possibility  of  plastic  strains  is  recognized,  four  questions 
ari  se: 


1.  Can  plastic  strains  occur? 

2.  If  they  occur,  what  will  be  their  relative  algebraic  values? 


^  l ..III* I-  •  <,  « 


3.  If  they  occur,  what  v/ill  be  their  actual  algebraic  values? 

4.  Will  they  occur? 

Obviously,  Question  2  is  a  subset  of  Question  3.  The  reasons  for  listing 
the  two  questions  separately  are  the  mathematical  and  physical  conditions 
used  to  answer  them,  which  are  explained  below. 

The  mathematical  theory  of  elastoplasticity  presented  here  contains 
four  parts,  each  needed  to  answer  one  of  the  above  four  questions: 

1.  A  yield  criterion,  assumed  to  be  of  the  form 

fU.j,  eP.)  =  0  (D.5) 


satisfaction  of  which  is  a  necessary  condition  for  the  occurrence 
of  additional  plastic  strain  at  a  point. 

2.  A  plastic  potential  function,  g(ojj),  for  which 

def,  =  dA  ia.  (D.6) 

lj  ^-jj 


which  gives  the  relative  algebraic  values  of  the  plastic  strain 
increments,  i.e.  the  direction  of  the  plastic  strain  increment 
vector  in  stress  space.  Equation  (D.6)  is  called  a  flow  rule, 
and  the  scalar  constant  dx  is  determined  by  the  next  condition. 

3.  The  requirement  that  Equation  (D.5)  be  satisfied  not  only  at  the 
beginning  of  yielding,  but  throughout  yieldi ng  as  well,  so  that 


df  =  do  •  • 

3o..  1J 
1  J 


(D.7) 


Equation  (D.7)  is  called  the  consistency  condition,  and  yields 
the  value  of  dx  in  Equation  (D.6).  It  therefore  permits 
calculation  of  the  actual  algebraic  values  of  the  plastic  strain 
i ncrements. 

4.  The  requirement  that  the  calculated  plastic  strain  increments 
lead  to  a  positive  plastic  work  increment, 

dWp  =  o-jjdefj  >  0  (D.8) 


139 


Equation  (D.8)  is  called  the  dissipation  condition.  If  it  is  not 
satisfied,  then  additional  plastic  strain  does  not  occur  at  a 
point,  in  which  case  all  strain  increments  are  elastic. 

Equations  (D.8),  (D.6),  (D.7),  and  (D.8)  have  been  written  assuming 
one  yield  criterion  (or  yield  surface),  and  one  plastic  potential 
function.  There  can,  however,  be  more  than  one  yield  surface,  and  an 
equal  number  of  correspondi ng  plastic  potential  functions.  If  this 
happens,  then  the  above  four  equations  apply  to  each  active  yield 
surface.  Thus,  if  m  yield  surfaces  are  active,  there  will  he  a  set  of 
plastic  strain  increments  for  each  active  yield  surface,  the  values  of 
which  are  determined  by  4m  equations  (counting  Equation  (D.6)  as  one 
tensor  equation. ) 

The  stress  tensor  a.,  contains  nine  elements, 

0 


but  only  six  are  independent  because 

aji  =  °ij  (D-10 

Each  stress  component,  <7^.,  can  be  expresed  as  a  function  of  the 
three  principal  stresses,  a^,  and  the  nine  direction  cosines 

of  the  three  principal  stress  axes  with  respect  to  the  arbitrary  Cartesian 
axes  used  to  define  the  <7^.  However,  if  a  unit  vector  pointing  in  the 
direction  of  the  ith  principal  stress  axis  is  e . ,  then  because  the  three 
principal  stress  axes  are  orthogonal,  we  have 


(D.ll) 


I4n 


Equation  (0.11)  represents  six  independent  scalar  equations  involving 


the  nine  principal  stress  axis  direction  cosines.  Thus,  there  are  only 
three  independent  principal  stress  axis  direction  cosines.  Let  them  be 
■*1 ,  and  .  Then  we  can  write 

’ij=  TijUl*^*"3;  al*a2*a3)  (D 


If  a  material  is  isotropic,  the  dependencies  of  the  yield  function, 
f,  in  Equation  (D.S)  and  of  the  plastic  potential  function,  g,  in 
Equation  (0.6)  on  the  principal  stresses  ,  <7,,  and  ace 
independent  of  the  orientation  of  the  principal  stress  axes  with  respect 
to  the  coordinate  axes.  This  means  that  not  only  do  Oj,  a,,,  and 
not  appear  in  the  expressions  for  f  or  g,  but  also  the  stress  functions 
wk i c h  do  appear  in  those  expressions  are  insensitive  to  subscript 
interchanges,  i.e.,  they  are  symmetric  functions  of  oj,  and  ny 
The  total  stress  invariants  Ij,  I,,  and  I3  satisfy  the  required 
conditions  of  symmetry.  They  are,  from  Appendix  A: 
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(A. 29 


(A.  30 


Equation  (0.6)  gives  the  relative  values  of  the  plastic  strain 
increments,  from  which  the  relative  values  of  the  principal  plastic  strain 


increments  and  the  orientation  of  the  principal  plastic  strain  increment 
axes  can  he  determined.  Since  the  plastic  potential  function,  g,  is  a 


function  of  the  three  total  stress  invariants,  I ^ ,  I^,  and  I ^ »  given 
by  Equations  (A. 28),  (A. 28),  and  (A. 30),  we  have 

g  =  gdj.Ip.I3)  (D.  13) 

so  that  Equation  (D.6)  can  he  written  in  the  form 
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Now  Equations  {A. 28)  and  (A. 29)  yield 
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To  obtain  the  derivatives  al^/ao. .  we  first  expand  the  expression  for 
I3  in  Equation  (A. 30).  Referring  to  the  development  of  Equation  (A. 27), 
we  obtain 

(D. 17) 
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Differentiation  of  Equation  (D.17)  now  yields 


w^ere 
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Z  -j  i  =  cofactor  or  signed  minor  of  a-jj,  arising  in  the  Laplace 
expansion  of  j^|  ,  which  is 

N  =  r3  =  ?  7ij  ):ij  =  7  Tr  $ 


(0.19) 


The  compact  form  of  Equations  (D.18)  is 
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Substitution  of  Equations  (D.15),  (D.16),  and  (D.20)  into  Equation  (D.14) 
yields 


(a..  -  1,5.  .)  +  iS-  v.  .1 


llj 


al,  "i J 


(D.  21 


Since  z  bas  the  same  principal  axes  as  does  n_,  it  follows  from 
Equation  (D.21)  that  de_p  also  has  the  same  principal  axes  as  does  n_. 

This  condition  is  a  consequence  of  the  assumption  of  material  isotropy, 
and  not  an  independent  assumption. 

A  convenient  assumption  concerning  the  dependence  of  the  yield 
function,  f,  in  Equation  (D.6)  on  plastic  strain  is  that  f  is  a  function 
of  stress,  a..,  and  plastic  work,  Wp,  where  plastic  work  is  in  turn  a 

'  J 

function  of  plastic  strain  [Malvern  (1960:367)1. 
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The  plastic  work  increment  in  Equation  (D.8)  can  be  written  in  the  form 
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so  that 
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1 1 


Substitution  of  Equation  (0.21)  into  Equation  (D.23)  yields 
d«P  .  c,  .dtP.  .  d>  [if.  Trio)  .  if-  [irl/l  -  if 


where 


=  dx 


ifr  'i  *  2  ifr  -2  *  3  If-  I. 


=  hdx 


+  fr^  h  Tr(^  - 


(D.25) 


h  =  I.  ♦  2  if-  I~  +  3  IS-  K 
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The  above  expression  for  h  agrees  with  Euler's  theorem,  since  the 
invariants  I^,  I^,  and  I3  are  homogeneous  functions  of  degree  1,  2, 
and  3,  respectively  LSokolnikoff  and  Redheffer  (1966:325)]. 

Equation  (0.13),  which  assumes  the  plastic  potential  function  depends 
explicitly  on  the  total  stress  invariants,  1^,  i^,  and  I3,  gives 
reasonably  good  results  for  stress  or  strain  paths  involving  mainly 
volumetric  compression.  However,  shear  strength  data,  plotted  in  the 
octahedral  plane  as  described  in  Appendix  C,  suggest  that  a  plastic 
potential  function  for  stress  or  strain  paths  involving  significant  shear 
deformation  is  better  assumed  to  depend  explicitly  on  the  first  total 
stress  invariant,  1^,  and  the  second  and  third  deviator  stress 
invariants,  J2,  and  J3,  so  that 
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Equation  (D.6)  then  yields 


deP .  =  dx  =  dx 

1J  3aij 


ag  3ll  +  _ag  SJ2  3skl  +  _ay_  aJ3  3Sk1 
3I1  3<Jij  3vJ 2  3skl  30ij  aJ3  3Skl  3oij 


(0.28) 


/'I 


145 


and  by  analogy  with  Equations  (D.16)  and  (D.2u) 
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and  since  Equation  (A. 34)  can  be  written  in  the  form 
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we  have 
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Substitution  of  Equations  (0.15),  (D.29),  (D.30),  and  ( D. 32 )  into 


Equation  (D.28)  now  yields 


The  plastic  work  increment  in  Equations  (D.8)  or  ( D . 2 3 )  is  therefore 


wh  e  re 
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Comparison  of  Equations  (D.26)  and  (D.35)  shows  that 


2  l2  +  3  al3  l3  ~  2  aJ2  J2  +  3  "iJj  J3 


(D.34) 


(D.35) 


( D.36) 


Finally,  the  plastic  work  increment,  dk/P,  can  be  expressed  as  the  sum 
of  a  volumetric  term  and  a  deviatoric  or  distortional  term.  This  is  done  by 
expressing  the  stress  components,  oij,  and  plastic  strain  increments, 
del-j,  as  the  sums  of  their  volumetric  and  deviatoric  components. 
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The  expression  for  the  plastic  work  increment  can  now  be  written  in  the  form 
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The  first  term  in  Equation  (D.39)  is  the  volumetric  plastic  work  incre¬ 
ment;  the  second  term  is  the  deviatoric  or  distortional  plastic  work  increment 
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VECTOR  REPRESENTATION  OF  A  GENERAL 
STRESS  OR  STRAIN  STATE 


Consider  the  two-dimensi onal  rectangular  Cartesian  coordinate  system 
(X,Y)  shown  in  Figure  (E.la).  The  X  and  Y  axes  define  a  two-dimensional 
vector  space,  and  the  line  x  =  y  is  a  one-dimensional  vector  subspace. 
Point  P  in  the  one-dimensional  subspace  can  be  located  by  the  single 
coordi  nate  s  =  ^2  x . 

In  three  dimensions  a  similar  situation  exists,  as  shown  in 
Figure  (E.lb).  The  X,Y,  and  Z  axes  define  a  three-dimensional  vector 
space,  and  the  plane  x  =  y  containing  the  Z  axis  is  a  two-dimensional 
vector  subspace.  Point  P  in  the  two-dimensional  subspace  can  be  located 
by  the  two  coordinates  (s  =  Jz  x.,  z). 

Generalization  of  the  above  analysis  for  an  n-dimensional  vector 
space  is  straightforward.  If  an  n-dimensional  vector  space  is  defined  by 
the  coordinates  x^  (i  =  l,n),  then  the  hyperplane  xi  =  Xj  (i  ±  j)  is 
an  (n-l)-dimensional  vector  subspace.  A  point  in  that  subspace  can  be 
located  by  the  coordinates 
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which  has  nine  components.  The  nine  stress  components  can  be  used  to 
define  a  nine-dimensional  vector  space,  in  which  the  location  of  a  point 
is  specified  by  a  vector  having  nine  components. 


{°9}  = 


(E.2) 


However,  only  six  stress  components  in  Equations  (E.l)  and  (E.2)  are 
independent  because  the  stress  matrix  is  symmetric  according  to 
Equation  (A. 20) ,  i .e. , 

cji  =  aij 

so  that 


°21  =  a12 
a32  =  a23 
al 3  =  a31 


(E.4a) 

(E.4b) 

(E.4c) 


Therefore,  only  the  six-dimensional  vector  subspace  defined  by 
Equations  (E.4),  of  the  nine-dimensional  vector  space  defined  by 
Equation  (E.2)  will  ever  be  needed.  Thus,  a  general  stress  state  can  be 
specified  by  the  six-dimensional  vector 
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Note  carefully  that  c^,  and  in  Equation  (E.5)  are  not 
necessarily  principal  stresses,  but  simply  the  first  three  elements  of  the 
column  vector  jaj  ,  which  represent  the  normal  stresses  a^,  a^,  dnd  °33* 
Expressions  for  the  total  and  deviator  stress  invariants  can  be 
developed  using  the  notation  of  Equation  (E.5).  The  total  stress 
invariants,  I,,  1 2 »  and  ^3*  are 

*1  =  °1  +  °2  +  CT3  (E’6) 
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I2  =“(cr^a2  +  ct2°3  +  a3°l^  +  2  ^°4  +  °5  +  °6  ^ 


t  *  1  1,2.  2.  2. 

3=al°2°3  O405a6  ‘  1  °1°5  °2  °6  °3  a4 
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then  the  deviator  stress  vector  can  be  written  in  the  form 


and  the  derivatives  of  the  deviator  stress  invariants  are 


(E.l 


(E.2 


where 
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■jj  {m}  {m}  T )  (S)  =  { S> 


In  a  similar  manner,  a  general  strain  state  can  be  represented  by  a 


six-dimensional  column  vector.  The  total  strain  vector  is 


APPENDIX  F 


INCREMENTAL  FLEXIBILITY  MATRIX  FOR 
STRESS  CONTROL 


Consider  a  rate  independent  elastoplastic  model  having  two 
independent,  strain  hardening  yield  surfaces  and  two  corresponding  plastic 
potential  functions.  Using  nomenclature  similar  to  Lade's,  one  yield 
surface  will  be  called  the  compressive  yield  surface  and  the  other  the 
expansive  yiel d  surface.  Four  incremental  deformation  modes  are  possible 
for  such  a  model : 

1.  Both  yield  surfaces  active,  which  implies  both  compressive  (C) 
and  expansive  (P)  yielding,  plus  elastic  (E)  response  (the  ECP 
mode) . 

2.  Only  the  compressive  yield  surface  active,  which  implies 
compressive  yielding  plus  elastic  response  (the  EC  mode). 

3.  Only  the  expansive  yield  surface  active,  which  implies  expansive 
yielding  plus  elastic  response  (the  EP  mode). 

4.  Only  elastic  response  (the  E  mode). 

For  such  a  model  a  plastic  strain  increment  has  two  components: 
compressive  and  expansive. 

.  p  .  pc  ,  .  pp 

deij  *  dTj  *  d‘fj 

Thus,  a  total  strain  increment  is  the  sum  of  three  components: 
elastic,  compressive  plastic,  and  expansive  plastic. 


(F.l) 


de . .  =  dEe.  +  deP.  =  dEP.  +  dcPP  +  dEPP 

ij  ij  ij  ij  ij  ij 

The  compressive  and  expansive  plastic  work  functions  are 
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Equation  (F.9)  into  Equation  (F.7),  and  Equation  (F.10)  into 
Equation  (F.8)  yields 


Equations  (F.13)  and  (F.14)  can  be  written  together  in  the  form 
•df  1  =  FT  {da)  -  rDj  {dA }  =  {0} 


where 


where  the  el astopl astic  incremental  flexibility  matrix,  HleP,  is  given  by 
the  expression 


Hep  =  He  +  G  rDj  _1  FT 


(F.25) 


Now  Equations  (F.13)  and  (F.14)  are  uncoupled  when  {da)  is  known,  so  that 


(F.  26) 


(F. 27) 


which  is  the  sane  result  obtained  by  expanding  Equation  (F.19). 
Equation  (F.25)  can  thus  be  written  in  the  form 


(F.28) 


When  the  EC  mode  is  active,  d>p  is  set  equal  to  zero  and  dxc  is 
calculated  by  Equation  (F.2fi),  as  before.  When  that  result  is  positive  the 
elastoplastic  incremental  flexibility  matrix  is  obtained  by  deleting  the 


expansive  tem  from  Equation  (F.28),  which  yields 


Hep  =  H ' 


T 


(F.  29) 


When  the  EP  mode  is  active,  dxc  is  set  equal  to  zero  and  dxp  is 
calculated  by  Equation  (F.27)  as  before.  When  that  result  is  positive  the 
elastoplastic  incremental  flexibility  matrix  is  obtained  by  deleting  the 
compressive  term  from  Equation  (F.28),  which  yields 


When  the  E  mode  is  active,  both  dx  and  dx„  are  set  equal  to  zero, 
and  the  elastic  incremental  flexibility  matrix  is  obtained  by  deleting  both 
the  compressive  and  expansive  terms  from  Equation  (F.28),  which  yields 

Hep  =  Hg  (F.31) 

The  logic  for  deciding  which  incremental  deformation  mode  is  active, 
under  stress  control,  is  discussed  in  Appendix  H. 


APPENDIX  G 


INCREMENTAL  STIFFNESS  MATRIX  FOR 
STRAIN  CONTROL 


The  material  in  this  appendix  builds  on  the  development  in  Appendix  F 
When  the  ECP  mode  is  active,  the  consistency  equations  for  a  two 


yield  surface  model  take  the  form  of  Equation  (F.15). 

{df }  =  FT  {do}  -  rDj  {dx}  =  {0} 

However,  under  strain  control  what  is  known  is  not  the  stress 


(F.15 


increment,  {do},  but  the  total  strain  increment,  {de } .  Thus,  in  place  of 


Equation  (F.21)  we  write 


{do}  =  C  e {d  e e } 


(G.l) 


where 


£e  =  elastic  incremental  stiffness  matrix,  defined  by 
Equation  (J.32). 

In  addition,  we  write  Equation  (F.20)  in  the  form 

(deej  =  {de}  -  {deP} 


(G.2) 


and  again  use  the  flow  rule  to  obtain  an  expression  for  the  plastic  strain 


i ncrement. 


{dep}  =  G  {dx} 

Substitution  of  Equations  (G.l),  (G.2),  and  (F.22)  into 
Equation  (F.15)  yields 

{df}  =  FTCe  jj de}  -  G  {dx}]  -  rDj  {dx}  =  {0} 

or  r  r-  T 

(F  CeG  +  'Dj )  { dx }  =  F  Ce  {de} 


(F.22 


(G.3) 


(G.4) 


Solving  Equation  (G.4)  for  { dx>  yields 

(dx)  =  (FTCeG  +  rDj)'1  FTCe  f  de }  (G.5 

Provided  dxc  and  dXp  are  both  positive,  we  substitute 
Equations  (G.2),  (F.22)  and  (G.5)  into  Equation  (G.l)  to  obtain  an 
expression  for  the  stress  increment,  {dn},  in  terms  of  the  prescribed 
total  strain  increment,  {del. 

{da}  =  Ce{dEe}  =  Ce[(de}  -  G  {dx}] 

=  Ce[fde}  -  G  (FTCeG  +  rDj)-1  FTCe{d£}] 

=  -  CeG  (FTCeG  +  rDJ)"1  FTCe]fdel 

=  Ccp  {de}  (G. 6 

where  the  elastoplastic  incremental  stiffness  matrix  for  the  ECP  mode  is 
defined  by  the  expression 

cep  =  Ce  -  CeG  (FTCeG  +  IDJ)"1  FTCe  (G.7 

The  fact  that  _Hep  as  defined  by  Equation  (F.25)  and  £ep  as 
defined  by  Equation  (G.7)  are  the  inverse  of  one  another  can  be  proven  by 
showing  that 

HepCep  =  CepHep  =  (G.8 

In  verifying  Equation  (G.8)  use  is  made  of  the  fact  that 

HeCe  =  CeHe  =  l  (J  .3 

When  the  EC  mode  is  active,  dx^  is  set  equal  to  zero  so  that 


Equation  (G.3)  reduces  to 


Solving  Equation  (G.10)  for  dxc  yields 


(G. 


Provided  d*c  is  positive,  we  substitute  Equations  (G.2),  (F.22), 
and  (G.ll)  into  Equation  (G.l)  to  obtain  an  expression  for  the  stress 
increment,  (do},  in  terms  of  the  prescribed  total  strain  increment,  { de } 

t  d<i }  = 


i 


=  £ep  {  de)  (G. 

where  the  elastoplastic  incremental  stiffness  matrix  for  the  EC  mode  is 
defined  by  the  expression 


The  fact  that  as  defined  by  Equation  (F.?h)  and  £ep  as 
defined  by  Equation  (G.13)  are  the  inverse  of  one  another  can  be  proven  by 
showing  that 

HepCep  =  CopHep  =  I  (G.8) 


again  using  Equation  (J.31). 

When  the  EP  node  is  active,  dxc  is  set  equal  to  zero  so  that 
Equation  (G.3)  reduces  to 
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_£  ce  — ^  + 

,3o  J  —  L  3a  J 


°22  dAp  =  l 3a 


Ce  {de } 


(G.  15 


Solving  Equation  (G.15)  for  dxp  yields 
r.na  T 


dA  = 

P 


bn  ^{de} 

r  ,n  T  f  a 

Pfn(  e  J3gol 

i — h  c  1 — £  +  I 

L3a  J  —  L3a  J 


(G.  16 


Provided  dXp  is  positive,  we  substitute  Equations  (G.2),  (F.22), 
and  (G.16)  into  Equation  (G.l)  to  obtain  an  expression  for  the  stress 
increment, (dc;,  in  terms  of  the  prescribed  total  strain  increment,  {de} . 

!  da}  =  Ce  {de6}  =  c6  [{de}  -  dx_] 


=  CGp{dc'(  (G.  17 ) 

where  the  elastoplastic  incremental  stiffness  matrix  for  the  EP  mode  is 
defined  by  the  expression 


The  fact  that  Hep  as  defined  by  Equation  (F.30)  and  Cep  as 
defined  by  Equation  (G.18)  are  the  inverse  of  one  another  can  be  proven  by 
showing  that 

HepCep  =  CepHep  =2  (G.8) 

again,  using  Equation  (J.31). 

When  the  E  mode  is  active,  both  dx_  and  dx_  are  set  equal  to 

L  p 

zero,  so  that  Equation  (G.2)  reduces  to 

{dEe}  =  {dE}  (G. 19) 

and  therefore  Equation  (G.l)  takes  the  form 


{ da)  =  CefdE)  =  Cep{dE}  (G.20) 

where  the  elastoplastic  incremental  stiffness  matrix  for  the  E  mode  is 
simply  the  elastic  incremental  stiffness  matrix. 

Cep  =  Ce  (G. 21 ) 
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APPENDIX  H 


INCREMENTAL  DEFORMATION  MODE  LOGIC  FOR  STRESS  CONTROL 


The  material  in  this  appendix  builds  on  the  development  in  Appendix  F. 
Equations  (F.26)  and  (F.27)  can  be  written  in  the  form 


dfc 

dxc  =  13 — 
c  un 


df ' 

dx  = 

P 


(H.  1 ) 


(H.  2) 


where 


df' 


(H.  3) 


d,P 


T 

fd«l 


( H .  4 ) 


When  D^  iS  positive,  implying  strain  hardening  behavior  of  the 

l 

compressive  yield  surface,  dfc  must  be  positive  for  dxc  to  be 

I 

positive.  The  same  argument  applies  to  D2 ?,  dfp,  and  dxp. 

Therefore,  in  view  of  Equations  (F.9),  (F.10),  (F.ll),  and  (F.12),  and  the 
dissipation  condition  as  stated  in  Equation  (D.8),  the  incremental 
deformation  mode  logic  for  a  strain  hardening  material  under  stress 
control  is  as  shown  in  Table  (H.l). 

When  both  the  compressive  and  expansive  yield  criteria  are  satisfied, 
i.e.,  when  the  stress  point  lies  at  the  intersection  of  the  compressive 
and  expansive  yield  surfaces  (the  corner  conditions  fp  =  fp  =  0)  the  above 
incremental  deformation  mode  logic  can  also  be  represented  graphically,  as 


shown  in  Figure  (H.l).  The  figure  is  drawn  in  .he  plane  of  principal 

stress  space  containing  the  yield  surface  gradient  vectors  Vf^  and  Vfp, 

because  Equations  (H.l),  (H.2),  (H.3),  and  (H.4)  show  that  the  quantities 

dx  and  dx  are  proportional  to  the  dot  products  of  da  with  "vf'  and  vf ' , 
c  p  c  p 

respectively.  It  follows  that  dx  and  dx  are  determined  by  the  component 

c  p 

of  do  in  the  plane  containing  "vf^  and  vf^.  Thus  if  we  set 


vf ;  x  vf' 

C  p 


(H.  5) 


then  n  will  be  a  unit  vector  perpendicular  to  the  plane  containing  vf^  and 

vf ' .  And  i f  we  set 
P 

_ ★  _  _  _ 

da  =  da  -  (da  •  n)n  (H.6) 


then  da*  will  be  the  component  of  "Ha  in  the  plane  containing  vf^  and 

vf'.  It  is  da*  which  determines  dx„  and  dx„  for  the  corner  condition, 

P  c  p 

as  shown  in  Figure  (H.l).  The  component  of  do  normal  to  the  plane 

containing  "vf'  and  vf',  i.e.,  ( a"  •  7T)1T,  causes  only  elastic  deformation, 
c  p 

The  equations  of  Appendix  F  can  now  be  generalized  for  the  corner 
condition  by  using  the  ramp  function,  R(x),  defined  in  Figure  (H.2)  as 


R(x) 


<5  { 2 )  dz  dy  =  x 
=  0 


(x  >  0) 
(x  <  0) 


(H.7) 


so  that  one  incremental  flexibility  equation  yields  the  total  strain 
increment  for  any  stress  increment  at  the  corner.  The  generalization  of 
Equation  (F.24)  is 


f3gcl 

1  1  ] 

P9d1 

i.3a  J 

>  R(df')  + 

1  c  °22  ! 

(H.  8) 


Equation  (H.8)  applies  for  all  four  incremental  deformation  modes  at 
the  corner. 

Negative  values  of  or  D22,  implying  strain  softening,  result 
in  lack  of  uniqueness  under  stress  control.  For  example,  if  is 

I 

negative  a  stress  change  for  which  df  is  negative  can  be  caused 
either  by  expansive  yielding  or  elastic  unloading,  and  there  is  no  way  to 
distinguish  between  the  two  under  stress  control. 

A  computational  problem  which  needs  to  be  addressed  is  how  to  avoid 
violation  of  a  yield  surface  which  is  inactive  at  the  start  of  a  stress 
increment.  If  the  stress  point  lies  beneath  a  yield  surface  at  the  start 
of  a  stress  increment,  what  assurance  is  there  that  the  stress  increment 
will  not  be  so  large  for  the  distance  to  the  inactive  yield  surface  so 
small)  that  the  stress  increment  "punches  through"  the  inactive  yield 
surface  The  answer  is  "none"  unless  a  restriction  on  the  magnitude  of 
the  stress  increment  vector  is  established  to  prevent  yield  surface 
violation  in  the  EC,  EP,  and  E  modes.  When  the  compressive  yield  surface 
is  inactive,  the  restriction  on  (da)  is 


F'  =  [— I  T  (dal  <  f"  -  f'  =  -  f 
C  Ua  J  —  C  C  C 


df;=l7fj  {da]  -  c  fc  =  '  fc  (H-9) 

and  when  the  expansive  yield  surface  is  inactive,  the  restriction  on  (dal  is 

Ki  T 

df  *  =  — E  {dal  <  f  -  f’  =  -  f  (H.  10) 


Table  ( H . 2 )  shows  that  for  each  combination  of  initial  conditions  and 


incremental  deformation  node  a  distinct  set  of  four  conditions  must  be 
satisfied  to  ensure  positive  energy  dissipation  and  prevent  yield  surface 
violation.  Each  set  is  a  subset  of  the  following  four  general  conditions: 

df  <  -  fc  (H. 11) 


df  <  -  f 
P  -  P 


dx„  >  0 

c  — 


(H. 12) 
(H. 13) 


dx  >0 
P  - 


(H. 14) 


Expressions  (H.ll)  and  (H.12)  ensure  that  fr  and  f  are 
nonpositive,  and  Expressions  (H.13)  and  (H.14)  require  that  dxc  and  dxp 
be  nonnegative. 

The  restrictions  on  the  magnitude  of  the  stress  increment  vector 
defined  by  Equations  (H.9)  and  (H.10)  are  implemented  as  follows.  Assume 
the  expansive  yield  surface  to  be  inactive  and  the  compressive  yield 
surface  active  at  the  start  of  a  stress  increment,  i.e.,  f  =  0  and 
fp  <  0.  Then  if  the  ratio 


exceeds  1.0,  the  expansive  yield  surface  will  be  violated  unless  the  stress 
increment  is  reduced.  This  is  done  by  splitting  the  stress  increment  into 


The  stress  subincrement  fda)^,  is  just  sufficient  to  bring  the 
stress  point  into  contact  with  the  expansive  yield  surface.  The  remaining 
stress  subincrement,  (do^,  is  then  applied  assuming  the  expansive  yield 
surface  to  be  active  (ECP  or  EP  mode). 

A  similar  procedure  is  used  when  the  compressive  yield  surface  is 
inactive  and  the  expansive  yield  surface  active  at  the  start  of  a  stress 
increment,  i.e.,  f£  <  0  and  f  =  0,  using  the  ratio 


c 


When  both  yield  surfaces  are  inactive  at  the  start  of  a  stress 

increment,  i.e.,  f  <  0  and  f  <  0,  compute  both  C  and  c  Unless 
c  p  c  p 

both  ratios  are  less  than  1.0,  set 


{do}  =  {dolj  +  {da }?  (H>16 

where 

{  d°} i  =  \  {d°l  (H.20 

(da}2  =  (1-  I){da}  (H.21 

x,  -  larger  of  c  and  r  (H.22 

c  p 


The  stress  subincrement  {do}j,  will  bring  the  stress  point  into 
contact  with  one  yield  surface,  at  which  point  one  of  the  first  two  tests 


■v-'N'.'-v-'. 
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must  be  applied  to  see  whether  {da}  2  should  be  split  to  avoid  violation 
of  the  other  yield  surface. 

When  both  Cc  and  ^  are  greater  than  1.0,  e.g.,  when 
’■c  > 'p  >  1.0,  we  could  set 

'dal  =  —  fda)  ♦  - ~)  (dal  +  (  1  -  _L )  (do) 

c  p  c  p 

and  apply  each  stress  subi ncrement  in  succession  without  recomputing  the 
initial  conditions  at  the  start  of  the  second  and  third  stress 
subincrenents.  However,  this  is  not  recommended  because  it  is  numerically 
less  accurate  than  recomputing  a  new  set  of  initial  conditions  at  the  start 
of  each  new  stress  increment. 


TABLE  H.2 


RESTRICTIONS  IMPOSED  FOR  EACH  INITIAL  CONDITION 
AND  INCREMENTAL  DEFORMATION  MODE 


Initial 

Conditions 

Incremental 

Deformation  Mode 

fc 

fp 

ECP 

EC 

EP 

E 

0 

0 

dfc  =  0 

dfp  =  0 

QAC  >  0 

dAp  >  0 

dfc  =  0 

dfp  <  o 
dxc  >  0 
dAp  =  0 

dfc  <  0 

dfp  0 

dxc  =  0 
dxp  >  0 

dfc  1  0 

dfp£° 

dxc  =  0 
dxp  =  0 

0 

<  0 

n/a 

dfc  C 

dfp  1  -fP 
dxc  >  o 

dx  =  0 

P 

n/a 

d<  <0 

dfp  1  -fp 
dxc  =  0 

dxp  -  0 

1 

0 

n/a 

n/a 

dfc  <  'fc 

dfp  =  0 

dxc  =  0 

dAp  >  0 

dfc  <  -fc 

dfp  <  0 

dxc  =  0 

dAp  =  0 

<  0 

<  0 

n/a 

n/a 

n/a 

dfc  <  -fc 
afp  i  -fp 

P»c  .  0 

<hp  ,  0 

175 


&’Ai 


"Frv  -Ty  *  ■  V  > 


R(x) 


x  y 

R(x)  =  /  /  <5(z)  dzdy  =  x  (x  >  0) 

-OO  -CO 

=0  (x  <  0) 


Figure  H.2.  The  ramp  function 


APPENDIX  I 


INCREMENTAL  DEFORMATION  MODE  LOGIC  FOR  STRAIN  CONTROL 

The  incremental  deformation  mode  logic  for  a  strain  hardening,  two 
yield  surface  el astoplastic  model  under  stress  control,  as  presented  in 
Appendix  H,  is  simple  because  the  two  yield  criteria  are  defined  in  stress 
space.  For  the  corner  condition,  the  incremental  deformation  mode  is 
determinded  by  the  component  of  "da  in  the  plane  containing  Yf^  and  Yfp. 

All  stress  increment  vectors  are  possible,  and  each  stress  increment 
vector  leads  to  one  and  only  one  strain  increment  vector,  as  indicated  by 
Equation  (H.8). 

<"'>  -  f  (d’;  *  TS^  [It]  r  {IH  r  ,dV  ,H'8) 

However,  it  remains  to  be  seen  whether  there  are  some  strain  increment 
vectors  not  produced  by  any  stress  increment  vector,  or  whether  two 
different  stress  increment  vectors  can  produce  the  same  strain  increment 
vector. 

The  situation  under  strain  control  is  more  complicated  than  that 
under  stress  control,  because  each  incremental  deformation  mode  has  a 
different  incremental  stiffness  matrix.  However,  the  conditions  listed  in 
Table  (H.2)  also  apply  under  strain  control,  and  what  must  be  done  is  to 
express  those  conditions  in  terms  of  a  prescribed  strain  increment  rather 
than  a  prescribed  stress  increment. 

The  conditions  df£  =  o  and  dfp  =  0,  when  they  apply,  are  taken 
care  of  by  the  consistency  conditions,  Equations  (F.7)  and  (F.8).  The 


conditions  dx  =  0  and  dx  =  0,  when  they  apply,  are  taken  care  of 
c  p 

i 

arbitrarily.  It  is  the  conditions  dxc  >  0,  dxp  >  0,  dfc  <_  -fQ9 

i 

and  df  <_  -f  which  require  further  discussion. 

First,  consider  the  corner  conditions  f  =  f  =  0.  For  the  ECP 
’  c  p 

mode,  Equation  (G.5)  can  be  written  in  the  form 
{dx}  =  L'1  {df * } 


where 


fLi>!  (L2> 


FTC6G  +  rD 


f3fcl  T  e  [39cl 

tar]  £ebr]  +  Di 

K]  T  ce  K] 

L3a  J  —  Lao  J 


[!!ci  t  c.  m 

(.So  )  —  is  o  } 

faf'l  T  „  fagj] 

_E  Ce  -P  + 
Lao  J  -  Lao  J 


Djj  and  are  defined  by  Equations  (F.ll),  (F.12),  and  (F.17),  and 

I 

the  elastic  trial  matrix,  { dfe > ,  is  defined  as 


= 


FTCe  {de} 


The  two  columns  of  the  vertically  partitioned  l  matrix  in 
Equation  (1.2)  can  each  be  interpreted  as  a  two  component  vector,  and  can 
be  represented  graphically  in  the  plane  containing  both  of  them  (the 
L  plane),  as  shown  in  Figure  (1.1).  The  lengths  of  the  column  vectors, 

L,  and  Lo  are 


and  the  acute  or  obtuse  angle  between  them  can  be  found  by  computing  their 
cross  product. 


where 


L1  x  L2  " 


11 

L12 

■21 

L22 

0 

0 

=  Le- 


sin  i  e^ 


L  =  L 1 1L  22  '  L 2 1L 1 2 


L11 

L12 

L21 

L22 

f L11 L?2  '  L21 L12 )e3 


(1.6) 


(1.7) 


Therefore,  if  L  >  0  the  acute  or  obtuse  angle  from  to  L2  is 
nonzero  and  counterclockwise. 

Now,  consider  the  inverse  of  JL,  which  appears  in  Equation  (1.1),  and 
set 


r,  "u 

H12  i* 

r  Mi  1 

lJj 

i  M21 

M 

Ll_£jJ 

f  lL?2 _ 

'Ll? 

i 

1 

1  3 

_ 

LU, 

(1.8) 


The  relation  between  the  row  vectors  of  the  horizontally  partitioned 
M  matrix  and  the  column  vectors  of  the  vertically  partitioned  L  matrix  is 


shown  in  Figure  (1.2).  In  order  for  M  to  be  the  inverse  of  L,  it  must  he 


that 


Mj  • 

Ll  =  1 

(1.9 

V 

L2  =  0 

(1.1 

V 

r!  =  ° 

(1.1 

V 

L0  =  1 

(1.1 

which  dictates  the  slopes  of  the  vectors  and  in  Figure  (1.2),  as 
well  as  their  lengths 


M 


1 

i =  r 


d.i 


(i.i 


Equation  (1.1)  can  now  be  expanded  to  yield 

dx  =  M.  •  cTf'  >  0 
c  1  e 


(I.i 


dx  =  M0  •  df *  >  0  (I • : 

p  2  e 

Equations  (1.15)  and  (1.15)  will  be  satisfied  if  and  only  if  the 
vector  dfg  lies  between  the  vectors  Lj  and  L?_  in  Figure  (1.2). 

That  zone  is  labelled  ECP  in  Figure  (1.3). 

When  df  is  parallel  to  L,,  Equation  (1.16)  gives 


A1  so,  in  that  case 

df '  df * 
e  ce 


•11 


so  that 


dx  =  M.  •  df'  =  M 
c  1  e 


rdf  df1 

ce  j-  _  ce 

L, ,  1  =  L. 


-11 


11 


Similarly,  when  dfg  is  parallel  to  L2,  we  have 


dxc  =  0 


df'  df' 
_ e  _  pe 


*22 


-  -  dfie 

dxn  =  M.  •  df'  =  |~ 

P  L.  0  L 


For  the  EC  mode  Equation  (G.ll)  can  be  written  in  the  form 


(1.18) 


(1.19) 


(1.20) 


(1.21) 


(1.22) 


r.>: 

V’  ■ 


df ' 


dx. 


ce 


>  0 


'll 


which  is  the  same  as  the  expression  obtained  when  dfp  -js  parallel  to 
Lj  in  the  ECP  mode.  When  >  0,  Equation  (1.19)  requires  that 


df '  >  0 

ce 


In  addition,  when  dxp  =  0,  Equations  (H.4),  (G.l),  (G.2),  and 


(F.22)  yield 


vv  vvT. 


dfpe  -  L21dxc 


dfpe  '  L21 


(1.24 


When  Ljj  >  o,  Equation  (1.24)  yields 

L11  dfi  •  L11  dfpe  -  L21  dfce  =  V  e  <  0 


(1.25 


The  zone  of  the  L  plane  in  which  both  Equations  (1.23)  and  (1.25)  are 
satisfied  is  the  EC  zone  in  Figure  (1.3). 

For  the  EP  mode,  Equation  (G.16)  can  be  written  in  the  form 


df 1 

dx  =  >  0 

p  L^T 

which  is  the  same  as  the  expression  obtained  when  df^  is  parallel  to 

L?  ln  The  ECP  mode.  When  >  0,  Equation  (1.22)  requires  that 

df'  >  0 
pe 

In  addition,  when  dxc  =  0,  Equations  (H.3),  (G.l),  (G.2),  and  (F.22) 
yield 

faf'*)  T  faf '  1  T  f  fag  *)  1 

df'  =  — ^  (da)  =  —  CeUde>  -  — ^  dx  > 

C  (,3a  J  C3a  J  —I  [do  J  pi 

=  df'  -  L19  dx 
ce  12  p 


(1.22 


(1.26 


df'  -  -r—  df'  <0 

ce  l_22  pe  - 


When  Lp2  >  0,  Equation  (1.27)  yields 


(1.27 


U,  df'  =  L00  df'  -  L,„  df'  =  M,  •  df’  <  0 


(1.28 


The  zone  of  the  L  plane  in  which  both  Equations  (1.26)  and  (1.28)  are 


satisfied  is  the  EP  zone  in  Figure  (1.3). 
For  the  E  mode, 

dx. 


{dA> 


1  dx 


=  {0} 


(1.29 


so  that  Equations  (H.3),  (H.4),  (G.l), 


(G. 2 ) ,  and  (F.22)  yield 

Ce  {  del  =  df^  £  0  (1.30 


and 

faf'l  T  faf'l  T 

dfP  =  Lioj  {do}  =  ^  {de}  =  dfpe  1  0  (1‘31 

The  zone  of  the  L  plane  in  which  both  Equations  (1.30)  and  (1.31)  are 
satisfied  is  the  E  zone  in  Figure  (1.3). 

Figure  (1.3),  which  is  called  a  polar  mode  check,  shows  both  the 
vectors  tj  and  L2  in  the  first  quadrant,  with  a  counterclockwise  acute 
angle  between  and  l_2.  However,  there  are  other  possibilities.  If 
we  assume  the  compressive  yield  surface  to  be  always  strain  hardening  so 
that  Ljj  >  0,  but  admit  the  possibility  of  strain  softening  for  the 
expansive  yield  surface  so  that  L22  might  be  negative,  then  there  are 
twelve  possible  relative  angular  positions  for  and  l2.  Because 
*“11  ?  Lj  must  iie  1n  The  first  or  fourth  quadrant.  But  L2  can 
lie  in  any  quadrant,  and  when  l?  lies  in  the  quadrant  opposite  (i.e., 
not  adjacent  to)  the  quadrant  containing  the  obtuse  angle  between 
Lj  and  L2  can  be  either  counterclockwise  or  clockwise.  The  twelve 
relative  angular  possibilities  for  Lj  and  l2,  together  with  their 


impact  on  both  uniqueness  and  completeness  of  the  incremental  deformation 
mode  solution  are  tabulated  below: 

QUADRANT  ANGLE 


ANGLE 

L}  to  L2 

UNIQUE 

COMPLETE 

CCW 

Yes 

Yes 

CCW 

Yes 

Yes 

cw 

Yes 

No 

cw 

No 

No 

cw 

No 

No 

CCW 

No 

Yes 

CCW 

Yes 

Yes 

CCW 

Yes 

Yes 

CCW 

No 

No 

cw 

No 

No 

cw 

Mo 

No 

cw 

No 

Yes 

In  the  above  table  a  unique  solution  is  one  having  no  overlap  between 
incremental  deformation  mode  zones,  and  a  complete  solution  is  one  for 
which  no  angular  zones  are  prohibited.  It  turns  out  that  the  only  cases 
for  which  the  incremental  deformation  mode  solution  is  both  unique  and 
complete  are  those  for  which 

Lu  >  0  (1.32 


l22  >  0 


L  =  LllL22  '  L12L21  >  0 


(1.33 

(1.34 


Now  consider  the  case  in  which  the  stress  point  lies  on  the 
compressive  yield  surface  but  below  the  expansive  yield  surface,  so  that 
fr  =  0  but  f  <  0.  Equation  (I. 19)  still  applies,  so  that 


There  is  no  danger  of  the  stress  point  punching  through  the 
compressive  yield  surface,  because  if  dxc  >  0  the  consistency  condition 
dfc  =  0  prevents  punch  through;  and  if  dxc  <_  0  it  is  because 

l 

dfce  _<  0,  so  that  the  stress  point  stays  on  or  pulls  away  from  the 
compressive  yield  surface.  However,  a  restriction  is  needed  to  avoid 
violation  of  the  expansive  yield  surface.  Equation  {I. 24)  can  be 
generalized  to  give 


df'  =  df' 


pe 


r—  R(df'  ) 


’ll 


ce 


(1.35) 


Then  if  the  ratio 


CP 


dfp 


(1.36) 


exceeds  1.0,  the  expansive  yield  surface  will  be  violated  unless  the 
strain  increment  is  reduced.  This  is  done  by  splitting  the  strain 
increment  into  two  parts  by  setting 

(d e  }  =  {dE}1  +  {de}2  (1.37) 


where 


{ d  c }  j  (dE } 

^  n 


(1.38) 


{de)2  =  (1  -  — )  (dt)  (1.39) 

CP 

The  strain  subincrement  {dc}^  is  just  sufficient  to  bring  the 
stress  point  into  contact  with  the  expansive  yield  surface.  The  remaining 
strain  subincrement,  {dE}2,  is  then  applied  assuming  the  expansive  yield 
surface  to  be  active  (ECP  or  EP  mode). 


eJI  - 
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A  similar  procedure  is  used  when  the  compressive  yield  surface  is 
inactive  and  the  expansive  yield  surface  active  at  the  start  of  a  strain 
increment,  i.e.,  f  <  0  and  f  =  0,  using  Equation  (1.22) 

L  p 

df ' 

dA  =  — >  0  (1.22) 

p  u22 


a  general ization  of  Equation  (1.27), 


dfc  •  df«  -  4f  R<dfpe> 


(1.40) 


and  the  ratio 


(1.41) 


When  both  yield  surfaces  are  inactive  at  the  start  of  a  strain 
increment,  i.e.,  fc  <  0  and  f  <  0,  compute  both 


(1.42) 


and 


C 

P 


(1.43) 


Unless  both  ratios  are  less  than  1.0,  set 

(del  =  { d  e  }  2  +  (de 


(1.44) 


where 


(de  },  =  —  {de } 

1  4 


(1.45) 
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fdc),  =  (1  -  -)  {del 
2  6 


(1.46) 


?  =  larger  of  6  and  6 


(1.47) 


The  strain  subincrement  (dc)^  w7]]  hrfng  the  stress  point  into 
contact  with  one  yield  surface,  at  which  point  one  of  the  first  two  tests 
must  be  applied  to  see  whether  (de  ^  should  be  split  to  avoid  violation 
of  the  other  yield  surface. 


I 


Incremental  deformation  mode  logic  for  a  strain 
hardening,  two  yield  surface  elastoplastic  model 
under  strain  control,  at  the  corner. 


APPENDIX  J 


ELASTIC  STRESS-STRAIN  EQUATIONS 


The  theory  of  linear  elasticity  assumes  the  stress  vector  {<?)  in 
Equation  (E.5)  and  the  strain  vector  { e T  in  Equation  (E.24)  are  related  by 
a  set  of  linear  equations  of  the  form 

•V-  -  hV  }  (J.l) 

where  the  elements  of  the  elastic  flexibility  matrix,  I4e,  are  constant. 

First  consider  the  work  done  when  two  successive  stress  increments, 
;d-  j  and  Ido)?.  are  added  to  an  existing  stress,  fal  The  increment 
of  work  must  be  independent  of  the  order  of  application  of  the  stress 
increments,  so  that 


Equation  (J.2)  yields 

fdr'j  {de'jj  =  (dol^  (deT^  =  {delj  {do},,  (J.3) 

and  from  Equation  (J.l)  we  have 


•;  dE'j  =  H 6  f  dolj 


(J.4) 


f  d e -2  =  H6  T  d c 1 ^ 


(J.5) 


Substitution  of  Equations  (J.4)  and  (J.5)  into  Equation  (J.3)  yields 


Substitution  of  Equations  (J.14)  and  (J.15)  into  the  left  hand  side 


With  H^e  defined  by  Equation  (J.23),  the  tensor  stress- strai n 
equations  represented  by  Equation  (J.l)  can  be  written  in  the  form 
1 


eij 


r  [(1  +  v)oij  '  ^kk6ij] 


(J.24) 


so  that  the  volumetric  strain  is 

1  -  2v 


-kk 


kk 


(J .25) 


and  therefore  the  deviator  strains  are 


Cij  EU 


ij  ■  ~r  6tj  =  r[(1  +  v)oij  •  vokk4ij]  ' 


1  -  2v  kk 


6ij 


(°ij  -  "7"  4ij)  =  ~ sij  (Jl26) 


Equations  (J.25)  and  (J.26)  can  be  inverted  in  the  form 


s .  .  = 


ij  ~  T  +”v  eij 


akk  =  1  -  2v  Ekk 


(J.27) 

(J.28) 


so  that 


c  +  °kk  ,  E 

O'*  -  S  .  .  *  — * —  6  •  •  — 

1 J  IJ  J  IJ 


t .  i  kk  ,  ^  t 

i£ij  -  -y  6ijj  3(1  -■■2v)  ekk6i j 


■  hi  *  (rhs)  ,J-29> 


which  means  that  the  elastic  stiffness  matrix,  C  ,  in  the  equations 


{0}  =  C  {c} 


(J.30) 


where 


» 


E(1  -  v) 


0  0 


0  0 


0  0 


(J.32) 


Equation  (J.31)  can  be  written  in  terns  of  two  quantities  fatnilar  in 
soil  mechanics,  viz. 


n  E(l-v) 

*  I  =  m', — vn — 


=  constrained  modulus 


(J .33) 


K  =  —  =  coefficient  of  lateral  stress  at  rest  (J.34) 

o  l-  v 


so  that 


1  Kn 

0  0 


0  0  0 


K  1 
o  o 


0  0  0 


Ce  =  M 


K  1 

0  0 


0  0  0 


1-K  0  0 

o 


(J .35) 


0  0  0 


0  l-Ko  0 


.0  0  0 


0  0  1-K. 


In  an  elastoplastic  model.  Equations  (J.l)  and  (J.29)  are  written  in 
incremental  form. 


{dee}  =  H6  (da) 


(J.36) 


(da)  =  Ce  (dee} 
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(J .37) 


r  r 


APPENDIX  K 


SPECIAL  EQUATIONS  FOR  THE  TRIAXIAL  TEST 


In  the  triaxial  test  the  normal  stresses  are  principal  stresses,  and 
the  lateral  principal  stresses  are  equal.  Consequently,  the  stress  matrix 
and  its  matrix  of  cofactors  reduces  to 


-7  0  0 

<1 

0  „  0 

r 

0  0 


(K.I) 


0  0 

a  0 


0  0 


If.?' 


and  the  total  stress  invariants  are 


h  =  qa  +  2nr 


1 2  =  ~qr'2r,i  +  ar^ 


L  o  ~  ^  > 

3  a  r 


y  y 


y.  .4; 


i  K  .  A  ) 


The  deviator  stress  matrix  and  its  matrix  of  cofactors  are 


h< ■  deviator  stress  invariants  are 


J 


i 


0 


(K 


(K 


(K 


+  o  . 
a  r 

— 1 


3  "2 


'  -  n 

a  r 


c  os 


J3 

(n a  '  nr\ 

.  ,  T 

\ 

, — 3  i 

J  J-  -  /J2VV2 

n  a  -  °r 

*3 

Vl) 

3 

3  =  sgn  (v  or) 


(K 


(K 


(K 


(K 


V,‘  .a.jsr.  of  axial  symmetry,  the  orientation  of  the  coordinate  axes, 
-  *r v.y  the  3x3  el astopol astic  incremental  stiffness  matrix,  CeP 


he  property  that 


The  3x3  el astopl astic  incremental  flexibility  matrix,  HeP)  ^as  the 
sane  properties. 


For  hydrostatic  compression 
aa  =  ar  =  11  OCT 

eKK 

Ga  ~  er  3 

and  the  3x3  incremental  stiffness  matrix  has  the  additional 

r  ep  _  fep  _  rep 
L11  ’  J22  '  l33 

reP  _  r<?P  _  rep  _  rep  _  rep  rep 
°12  =  13  “  21  “  ^23  “  °31  “  u32 

Therefore,  the  octahedral  normal  stress  increment  is 


(K.  19) 

(K.20) 


property  that 
(K.21) 

(K.22) 


'Cep  +  2Cep' 

d^0CT  =  daa  =  (C11  +  2C12)  dea  =  (  ~  )  deKK=  B'"PdeKK 


(K. 23) 


where  the  el astopl astic  incremental  bulk  modulus  for  pure  hydrostatic 
compression,  BeP,  is 


B 


CeP  +  2Cep 
ep  L11  ^l12 


(K.24) 


As  long  as  the  octahedral  normal  stress  increases  monotonical ly ,  the 
incremental  deformation  mode  is  the  EC  mode.  Unloading  and  reloading  are 
elastic. 

For  triaxial  compression  at  constant  cell  pressure,  starting  from 
isotropic  consolidation,  the  ECP  mode  is  active  as  long  as  (oa  -  ar) 
is  increasing.  Since 

do  =  0  (K.25) 

r 


DC 


I 


the  incremental  flexibility  equations  take  the  form 

d  rj 

de  =  dna  =  — 

a  11  a  ^-ep 


Hep 

de  =  H®P  d,  =  ^de  =  -v°Pdc 

r  21  a  uep  a  a 

H11 

where  the  el astopl astic  incremental  Young's  modulus,  Eep,  is 


rep  1 


and  the  elastoplastic  incremental  Poisson's  ratio,  vep ,  is 


(K  .2(1) 


(K . 27 ) 


(K. 28) 


(K.29) 


If  the  calculation  is  performed  under  axial  strain  control,  then 
Equation  (K.26)  is  written  in  the  form 


a  Hep 
M11 


=  Eep  de 


(K.30) 


Tne  above  equations  can  be  written  in  incremental  stiffness  form  as 


well,  by  setting 


d^a  =  ClPldEa+2C12der 


dv  =  0  =  c|P  dea  ♦  (Ce2P  ♦  C°2p)  der 


(K.31) 


(K .32) 


Cl 


- 


>*.y 


Equation  (K.32)  yields 


de  _ 


*ep 

'21 


^ep  +  pep 

l22  l23. 


de,  = 


-vep  de. 


where 


ep 


-op 

'21 


peP  +  rep 
l22  23 


Substitution  of  Equation  (K.33)  into  Equation  (K.31)  yields 


d"a  *  (cn  -  2^P  C!P)  d'a  ■  Eep  dea 


■where 


2C?p  Cep 
rep  _  rep  _  12  21 

'11  “  fep  fep 

l22  l23 


Conparison  of  Equations  (K.28)  and  (K.34),  and  (K.28)  and  (K.36) 


vep 


Hop  cep 

_21  _  U21 

Hep  "  fep  +  pep 
nll  l22  l23 


Eep 


1 


+  C 


ep 

23 


which  is  confirmed  by  forming  the  inverse  of  the  2x2  matrix  of 


(K.33) 

(K.34) 

(K.35) 

(K.36) 

shows  that 

(K.37) 

(K.38) 


coefficients  of  Equations  (K.31)  and  (K.32). 

For  a  general  triaxial  stress  path  test  in  which 


so  that 


da  -  da  . 
dq  _  _a _ r  _  1  -  a 

dp  "  da  +  do„  1  +  a 
a  r 


IK. 40 


the  incremental  deformation  mode  is  determined  by  the  methods  of 
Appendix  H,  since 

dfi-{£)T  ^ 


f3fnl  T 

dfn  =  {da) 

P  (. 3a  J 


where 


{da  }  =  da  <  a 


The  methods  of  Appendix  H  apply  no  matter  whether  the  calculations 
are  performed  under  axial  stress  or  strain  control,  because  the  stress 
path  is  prescribed  by  Equation  (K.39). 

The  incremental  flexibility  equations  take  the  form 

de  =  (H*!?  +  2a  H*P)  daa  (K 

a  11  I  l-  a 


dcr  .  [h^  .  .(Hg  * 


If  the  calculation  is  performed  under  axial  strain  control, 
Equations  (K .44)  and  (K.45)  are  written  in  the  form 


SE 


where 


8  = 


reP  rep 
C21  ~  a  C11 

CeP  +  reP  _  2n  Cep 
22  L23  L12 


Substitution  of  Equation  (K.52)  into  Equation  (K.4°)  yields 


(K. 53) 


K  -  <cn  *  26  cf2>  S 


If  the  calculation  is  performed  under  axial  stress  control, 
Equation  (K.54)  is  written  in  the  form 

dry. 


ie  ,  = 


^  reP  +  7a  reP 
L11  ^  L12 


For  constrained  (one-dimensional)  compression 


de„  =  0 


(K.54) 


(K . 55 ) 


(K . 56) 


A  constrained  compression  test  is  a  particular  strain  path  test,  so 
that  the  incremental  deformation  mode  is  determined  by  the  nethods  of 
Appendix  I.  The  incremental  stiffness  equations  take  the  form 

dr  =  CGp  de  =  MGp  dEn  (K.57) 

alia  a 


cep 

da  =  dP  de  =  da  =  KGP  da 
r  71  3  r^P  3  0  3 

where  the  incremental  el astopl astic  constrained  modulus,  Mep,  is 


(K. 50) 


Mep  -  Cep 
1  -  L11 


(K.59) 
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:v  ,:v,  r  S. 


•,.y 


and  the  incremental  elastoplastic  coefficient  of  lateral  stress,  is 


,ep  _21^ 
'o  rep 
L11 


( K . 60 ) 


For  constant  volume  compression 


dc^  +  2de  ^  =  0 


(K.6i: 


so  that 


der  =  -  2dea 


(K.62) 


A  constant  volume  test  is  also  a  particular  strain  path  test,  so  that  the 
incremental  deformation  mode  is  again  determined  by  the  methods  of 
Appendix  I.  The  incremental  stiffness  equations  take  the  form 


daa  "  (C1?  -  C1P2}  dea 


dar  "  [C21  -  1  (C2P2  *  CS>]d‘ 


(K.63) 


(K.64) 


so  that 


d n  a  a  da. 
r  a 


(K.65) 


where 


rep  1  ,rep  pepv 
U21  '  7  '22  '  V 


If  the  calculation  is  performed  under  axial  stress  control 


Equation  (K.63)  is  written  in  the  form 


de  =  - 2 - 

a  rep  rep 


For  a  general  triaxial  strain  path  test  in  which 


der  =  -  8  dea  (K.68) 

the  incremental  deformation  mode  is  determined  by  the  methods  of  Appendix  I 
The  incremental  stiffness  equations  take  the  form 

dca  =  (C*p  -  2  8  ceP)  de a  (K.69) 

dor  =  eg  -  8  (C®P  ♦  c|p)  de a  (K. 70) 

so  that 

da  =  a  da  C K . 7 1 ; 

r  a 

where 


C21  '  e(C22  +  C23} 

Cep  -  2s  Cep 
L11  L12 


(K.721 


If  the  calculation  is  performed  under  axial  stress  control. 


Equation  (K.69)  is  written  in  the  form 


APPENDIX  L 


TRANSIENT  RESPONSE  OF  A  THREE  ELEMENT 
VISCOELASTIC  MODEL 


Consider  the  three  element  viscoelastic  model  shown  in  Figure  (L.l) 
The  basic  equations  governing  the  model  response  are  [Bland  (1960:3)]. 

(L. 


Fi  =  F2  =  F 


X)  +  x2  =  X 


F1  =  klxl 


(L. 

(L. 


dx. 


F„  =  k„x„  +  iv 


2  “  2  2  2  dt 


(L. 


Substitution  of  Equations  (L.l),  (L.2),  and  (L.3)  into  Equation  (L.4) 
y  i  el  ds 


F  =  k. 


1  F! 

.  d 

rl 

r'2 

X 

1 

or 


k. 


r  ■,  '2  .  dx  n2  dF 

F  =  k2x  ‘  F  n2  It  *  Tc]  dt 


or 


dF  .  |  kl  *  k' 


dt 


Z,F.k  * 

Ho  /  ldt  n  n 


(L. 


Equation  (L.5)  can  be  written  in  the  form 


%  +  aF  =  g(t) 


(L  .6) 


Assuming  that 


x( 0- )  =  F(O-)  =  0 


(L. 


the  solution  to  Equation  (L.6)  is 


F(t)  =  f  g(y)e~di't~y'1  dy 


(L. 


Substitution  of  Equation  (L.8)  into  Equation  (L.10)  yields 


r  %. 


F(t)  =  k 


4*  + 


dY  +  —  f x(Y)e‘a(t'Y)dY 

l-o  Y  n?  i 


Ty  e 


(L. 


The  second  integral  on  the  RHS  of  Equation  (L.ll)  can  be  integrated 


by  parts  by  setting 


sc  that 


x|r)E-alt-'>  d,  ,  i  x<-,)Wa|t-’> 
a 


t  t 


CX  -d(t-v)  . 

37  e  d' 


0  0 


i  x[t)  .  /  OX  e-a(t-,)a 

a  /  fly 


( L . 16  ) 


ire  substitution  of  Equation  (L.16)  into  Equation  (L.ll)  yields 


1  k2 


I.?*’*11-'' 

K1 


dx 

d7  dY 


(L.17) 


If  the  Gi spl acement,  x(y),  is  a  Heaviside  step  pulse,  so  that  dx/dy 
is  a  Dirac-delta  function,  i.e. 


x( Y )  =  XH(y) 


(L .18) 


£  =  x  Sl> 
dy 


(L.19) 


where 


X  =  constant 


(1.20) 


(y  <  0) 
(y  =  0) 


(L  .21 ) 


(y  >  0) 


(L  .22 ) 


m 


i  *  Pi 

S-n-S 


/  / 


so  tl'.at  the  parameters  k  ^ ,  k^,  and  r,? 
kj  =  R(0+ ) 


R?  ( 0- ) 
R  (0+) 


can  be  calculated  from  the  equations 

( L .  30 

(L .  31 

(L  .32 


Note  that  Figure  (L.l)  shows  only  one  of  two  possible  viscoelastic 
nocel  s  containing  two  springs  and  one  cashpot  [Bland  (1960:3)],  and  that 
the  determination  of  the  parameter  values,  kp  k^,  and  is  conf i gurati on 
dependent  anc  therefore  not  unique. 
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APPENDIX  M 


YOUNG'S  MODULUS  FOR  A  HYPERBOLIC  STRESS-STRAIN  CURVE 


The  fact  that  a  hyperbolic  equation  of  the  form 


a  +  bx 


can  be  written  in  the  linear  form 


-  =  a  +  bx 

y 


(M.l) 


(M.2) 


has  long  been  used  to  obtain  an  empirical  formula  for  data  originating  at 
the  origin  and  rising  with  steadily  decreasing  slope  to  approach  a 
horizontal  asymptote  [Saxelby  (1913:138);  Running  (1917:39,53)]. 

Southwell  used  the  method  in  1932  to  determine  the  critical  load  of  an 
elastic  column  with  initial  curvature;  Gregory  used  it  in  1959  to  study 
structural  stability;  Kondner  and  Krizek  used  it  in  1962  to  fit  a 
hyperbola  to  footing  load--settlement  data;  Kondner  and  Zelasko  used  it  in 
1963  and  1964  to  fit  a  hyperbola  to  triaxial  compression  data  for  both 
sand  and  clay;  and  Duncan  and  Chang  used  it  in  1970,  also  to  fit  a 
hyperbola  to  triaxial  compression  data  for  soil,  as  a  means  of  obtaining  a 
general  expression  for  the  tangent  Young's  modulus  to  use  in  a  finite 
element  computer  code  [Southwell  (1932);  Southwell  (1936);  Gregory  (1959); 
Gregory  (1960);  Timoshenko  and  Gere  (1961);  Kondner  and  Krizek  (1962); 
Kondner  (1963);  Kondner  and  Zelasko  (1963);  Kondner  and  Zelasko  (1964); 
Duncan  and  Chang  (1970)]. 


i  i 


•  *  .  ^  '  .  ■ 


\  N  *.  *.  *. 
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For  the  case  of  cylindrical  compression  uncer  constant  confining 
(radial)  stress,  the  intermediate  and  minor  principal  stresses  are  equal 


and  constant,  so  that 


°2  =  a3  =  a3C 


Equation  (J.23)  therefore  yields 


dai  d<ai  '  °3C} 


■1  _  E, 


uhere 


Ey  =  tangent  Young's  modulus  for  a  nonlinear  stress- strai n  curve 
and,  therefore 


d(al  '  °3C ‘ 


If  the  relation  between  principal  stress  difference  and  axial  strain 
is  assumed  to  be  of  the  hyperbolic  form 


'1  '  3C 


then  it  foil ows  that 


e  ^  -*  0 


lim  (oj  -  °3C )  -  (oy  -  °3r^uit  =  F 


Thus,  the  parameter  a  in  Equation  (M.6)  is  the  reciprocal  of  the 
initial  slope  of  the  stress-strain  curve,  and  the  parameter  b  is  the 


reciprocal  of  the  upper  asymptotic  limit  of  the  principal  stress 
difference.  The  parameters  a  and  b  are  commonly  estimated  by  plotting 
E^/(o^  -  o-jq)  against  and  fitting  a  straight  line  to  that  plot, 
as  shown  in  Figure  (M.l),  because  Equation  (M.6)  has  the  linear  form 


°1  '  °3C 


=  a  +  be 


Note,  however,  that  the  values  of  the  parameters  a  and  b  obtained  by 
linear  regression  for  Equation  (11.9)  are  not  the  same  as  those  obtained  by 
nonlinear  regression  for  Equation  (M.6),  or  even  those  obtained  by  linear 


regression  for  the  equation 


°1  '  °3C 


=  a  i-  +  b 


(M.10 


because  the  error  functions  to  he  minimized  in  each  case  are  different. 
Note  also  that  if  we  set 


(al  -  °3C)ult  a 

- r. - =  3  =  e10 


( M .  1 1 


then  Equation  (M.6)  can  be  written  in  the  form 


Duncan  and  Chang  used  Janbu's  empirical  relation  between  E.  and 


'3C  [Janbu  (1963)],  which  is 


(M.  13) 


where 


p,  =  atmospheric  pressure 

a 

and  K  and  n  are  dimensionless  constants  often  obtained  by  fitting  a 
straight  line  to  a  plot  of 


versus  log.Q 


because  Equation  (M.13) 


has  the  linear  form 


lo9l0  li L\  =  lo9l0  K  +  n  1o910 

Va 


(M.14) 


Again  note,  however,  that  the  values  of  the  parameters  K  and  n 
obtained  by  linear  regression  for  Equation  (M.14)  are  not  the  same  as 
those  obtained  by  nonlinear  regression  for  Equation  (M.13). 

Duncan  and  Chang  also  used  the  Mohr-Coulomb  failure  criterion  to 
relate  measured  strength,  (a^  -  to  cell  pressure,  o^,  as  shown  in 

Figure  (M.2). 


,  s  2(c  cos  t  +  o nr  sin  <5) 

(01  ■  °3C  f  “  - . - r— £ - 

1  -  sin  t> 


(M.15) 


They  found  that  the  measured  strength,  (c^  -  °3c ) f »  tell  below  the 
fitted  asymptotic  strength,  (o^  -  a^)^,  ^  a  tactor  Rf,  so  that 

R* 


(a,  -  a  )  r  =  Rf  to,  -  O’)/-  ),,i  f 


(0.75  <  R,  <  1.00)  (M. 16) 


However,  the  above  observation  is  not  consistent  with  much  of 
Konclner's  published  data,  and  may  be  affected  by  the  hyperbolic  curve 


fitting  procedure  used. 

To  obtain  an  expression  for  Ey  as  a  function  of  -  o^,  w c 
write  Equation  (IT. 6)  in  the  form 

(ol  '  g3C} 

_ '  °3C  jul t 

£10  1  -  (ol  "  g3C; 

(ol  "  °3C)ult 


Differentiation  of  Equation  (M.17)  with  respc-ct  to 
yi el ds 

tie . 


°3C 


then 
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(M.l 


so  that 


d ( 0 1  '  c3C) 


de. 


Ej  .  Ei 


1  - 


( 0 1  "  °3C* 


1 2 


(or  03cJuit 


(M.l 


Substitution  of  Equations  (M.13),  (M.15),  and  (M.16)  into  Equation  (IT. 19) 
therefore  yiel ds 


E 


T 


°3C  \ 

Pa) 


Rf  (1  -  sin  0) (oj 
2  (c  cos  0  +  o~r 


(M.2 


What  Duncan  and  Chang  did  to  obtain  a  general  expression  for  Ey  for 
finite  element  code  use  was  assume  that  for  any  stress  path,  even  when 
o 3  is  not  constant,  the  relation  between  Ey,  Oj,  and  is  that 
obtained  by  replacing  in  Equation  (M.20)  by  o3,  to  obtain 


APPENDIX  N 

A  HYPERBOLIC  EXPRESSION  FOR  POISSON'S  RATIO 


In  a  triaxial  compression  test  at  constant  cell  pressure,  a  plot  of 
axial  strain  versus  radial  strain  for  stress  levels  well  below  the  failure 
level  might  appear  as  shown  in  Figure  (N.l)  [cf  Desai  and  Abel 
(1972:323)].  Such  a  plot  suggests  the  possibility  of  a  hyperbolic  fit,  of 


the  form 


which  means  that 


(N.l) 


(N.2) 


and  therefore 


del  (1  -  de2 ) 2  (1  -  d^)2 


where,  from  test  data  [Kulhawy,  Duncan  and  Seed  (1969)]  found  that 


vi  ’  f  “  G  -  F  log1(J 


The  parameter  d  is  obtained  as  the  slope  of  a  plot  of  -tg/tj 
versus  -e?,  since  Equation  (N.l)  has  the  linear  form 


-  f  ♦  <1  (-€,) 


The  parameters  F  and  G  are  obtained  from  a  semi  1 ogari thmic  plot  of 
“i  versus  1o9l0  (o3C/Pa) - 

Substitution  of  Equations  (M.17),  (M.13),  (M.1G),  and  (M.15)  into 
Equation  (N.3)  yields 


APPENDIX  0 


HYPERBOLIC  MODEL  FOR  CYCLIC  SIMPLE  SHEAR 


[Pyke  (1979:721)]  proposed  a  variation  of  the  hyperbolic  shear  model 
for  irregular  cyclic  loading,  in  which  the  simple  shear  stress-strain 
curve,  starting  at  the  latest  point  of  strain  reversal,  (■yc,tc),  is  a 
hyperbola  with  fixed  initial  slope,  and  fixed  upper  and  lower 

asymptotes,  and  -t  ,  as  shown  in  Figure  (0.1).  If  dy  >  0  after 
reversal,  hyperbola  PQ  applies,  having  the  equation 


T  -  T  = 
C 


gmax(y  -  Yc} 


MAX  , 

t  TT  (y  '  V 
y  c 


(dy  >  0) 


and  if  dy  <  0  after  reversal,  hyperbola  PR  applies,  having  the  equation 


T  -  T  = 
C 


gmax(yc  •  y) 


(dy  <  0) 


Equation  (0.2)  can  be  written  in  the  form 


-  O 


r-T^T-  <T  -  YC) 

y  c 


(dy  <  0) 


which  is  identical  to  Equation  (0.1)  except  for  the  sign  of  the  \ 


term.  Therefore,  Equations  (0.1)  and  (0.3)  can  be  combined  in  the  form 


APPENDIX  P 


YIELD  SURFACE  VIOLATION  CORRECTION  FOR  AN 
ELASTIC-PERFECTLY  PLASTIC  MODEL 

An  el astic-perfectly  plastic  model  does  not  strain  harden.  Instead 
the  yield  surface  remains  fixed,  and  so  it  is  often  called  a  failure 
surface.  The  failure  surface  is  concave  with  respect  to  the  origin  in 
stress  space,  and  its  intersection  with  any  octahedral  plane  is  a  closed 
curve,  with  its  center  on  the  hydrostatic  axis  and  having  six-fold  angular 
symmetry,  when  the  stress  point  lies  inside  the  failure  surface,  i.e., 
when  the  shortest  line  from  the  stress  point  to  the  hydrostatic  axis  does 
not  intersect  the  failure  surface,  the  material  responds  entirely 
elastically.  Plastic  strains  occur  only  when  the  stress  point  lies  on  the 
failure  surface;  and  the  stress  point  is  prohibited  from  going  beyond  the 
failure  surface  by  a  correction  procedure  which  will  now  be  explained. 

Consider  a  failure  surface  having  the  equation 


Its  octahedral  cross-section  is  a  circle  with  its  center  on  the 
hydrostatic  axis;  and  its  intersection  with  any  plane  containing  the 
hydrostatic  axis  is  a  scaled  plot  of  Equation  (P.l),  shown  in 
Figure  (P.l).  Assume  that  a  strain-controlled  elastic  trial  calculation 
starting  from  Point  1  in  Figure  (P.l)  has  moved  the  stress  point  to 
Point  3,  passing  through  the  failure  surface  at  Point  2.  The  objective  of 
the  correction  procedure  is  to  bring  the  stress  point  back  to  the  failure 
surface  at  constant  total  strain.  This  is  done  by  using  the  consistency 


condition  starting  from  Point  2,  and  a  (possibly  nonassociative)  flow  rule 
to  find  that  portion  of  the  total  strain  increment  which  is:  plastic, 
normal  to  the  plastic  potential  surface  at  Point  4,  and  just  sufficient  so 
that  the  remaining  elastic  strain  increment  causes  the  stress  increment  to 
terminate  in  the  failure  surface  (at  Point  5). 

The  consistency  condition  without  strain  hardening  reduces  from 
Equation  (D.7)  to  simply 

df  ■  (•£)  T  W  ■  °  (P-2) 


where 


and 


(P.3) 
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{m} 


3  f 
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II 


— ==.  (s)  [cf  Equation  (E.22)] 


(P.4) 

(E.14 

(P.5) 


(P.6) 


so  that  Equation  (P.3)  takes  the  form 


In  addition,  if  the  plastic  potential  function  is  of  the  same  form  as 
the  failure  criterion,  i.e., 

9  =  9  (I1  ’//'V  (p.6) 

and  i  f 


ag 

TT 
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=  9I 


(P.9) 


(P.10 


the  n 


(f)-, 


(p.n: 


how  if,  in  Equations  (J.27)  and  (J.28),  we  set 


3(T-2v') 


=  K 


=  2G 


(P.12) 

(P.13) 


where  K  and  G  are  the  elastic  bulk  and  shear  moduli,  respectively,  then 
Equation  (J.30)  can  be  written  in  the  form 


Ce  =  K  {m}  { m }T  +  2G  |l_  -  ^  {m}  { m 


=  (k  -  “)  {m}  {m  }T  +  2G  I 


(P.14) 


Then  with  a  view  toward  evaluating  Equation  (G.16),  Equations  (P.7)  and 
{ P . 1 4 )  yield 


can  be  calculatec  from 


Having  d*, 
the  equations 


the  adjustments  to  1^  anG 


CcKK  =  =  3 g j  <j a 


■,deh 


2  [j7 
'i  2 


(dlj)  =  -2K  dt^K  =  -  9Kgr  dx 


(P.2; 

(P.2- 


(P.2! 


{s}  •  20  {ce*3}  =  “GQjj  dx 


and  therefore  the  values  of  Ij,  /p2  ,  anc  { s  }  at  Poi nt  5  in 
Figure  (P.l)  are 


(Il)5  -  (I1>3  -  9KgI  dX 


(P.2i 


(P.2 


(P.2) 


(P.2' 


Note  that  the  derivatives  f j ,  fjj,  gj,  and  g^j  are  evaluated 
at  Point  4  and  not  Point  5.  Therefore,  the  values  of  1^  and 
calculated  by  Equations  (P.27)  and  (P.28)  may  still  define  a  stress  point 
lying  slightly  outside  the  yield  surface.  Therefore,  if 

f5  >  e  >  0  (P.3l 

where  e  is  a  specified  tolerance,  then  repeat  the  correction  procedure  by 
treating  Point  5  as  a  new  Point  3,  and  defining  a  new  Point  4. 


Figure  P.l.  Correction  procedure  for  elastic-perfectly 
plastic  failure  surface  violation. 


APPENDIX  0 

AFWL  ENGINEERING  MODEL  INCREMENTAL  PLASTIC  RESPONSE 


Figure  (Q.l)  shows  the  shear  failure  and  plastic  potential  surfaces, 
and  hysteretic  hydrostat  which  partially  define  a  modified  version  of  the 
AFWL  engineering  model  in  use  at  Applied  Research  Associates.  Each 
segment  of  the  failure  surface  has  an  equation  of  the  form 

f  ( 1 1  )  =  JjJ  -  (a  +  hlj)  =  0  (0.1 


and  the  von  Mises  plastic  potential  function  has  the  equation 

so  that  Equations  (P.4),  (P.5),  (P.9),  and  (P.10)  yield 

f!  =-b 

fII  "  1 
=  0 

9 II  -  1 


(0.2 


(Q.3 

(Q.4 

(0.5 

(0.6 


and,  therefore.  Equation  { P . 32 )  yields 


Substitution  of  Equation  (P.14)  into  Equation  (Q.7)  yields 


Now  the  total  strain  increment  can  be  written  in  the  form 


so  that  Equations  (Q.8)  and  (O.Q)  yield 


APPENDIX  R 


DRAINED  CAP  MODEL  AND  COMPUTATIONAL  ALGORITHM 


The  cap  model  is  an  isotropic,  rate  independent,  el astopl astic  model 
with  two  yield  surfaces.  The  shear  yield  surface  does  not  strain  harden, 
and  is  called  the  shear  failure  surface.  The  cap  yield  surface  both 
hardens  and  softens  in  response  to  plastic  volumetric  strain.  Both  yield 
surfaces  are  associative.  Figure  (R.l)  shows  the  two  cap  model  yield 
surfaces,  which  intersect  at  the  corner  point,  1^  =  k. 


The  shear  failure  surface  has  the  equation 


fjT  =  f(IJ 
1  2  1 


(Ij  >  T) 


(Ij  <  T) 


where  T  is  the  tension  cutoff,  and 


The  cap  yield  surface  has  the  equation 


(R.l) 


(R.2) 


(R.3) 
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(R.7) 
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(R.9) 

(R.10) 

(R.ll ) 

(R.12) 


When  k  <  0  a  von  flises  cap  transition  surface  is  defined  by  the  equation 

=  F(0,k)  =  f(k)  (k  <  Ij  <  0)  (R.13) 


and  therefore  when  <  0,  the  yield  criterion  is  the  lesser  of  the 
shear  failure  criterion  and  the  von  Mises  transition,  i.e., 

=  Fe  (Ij.k)  =  min  (fdj);  f(k j}  (T  <  0) 


(R.14) 


Within  the  yield  surfaces  the  material  is  hypoelastic  without 
hysteresis,  with  incremental  bulk  and  shear  unloading/rel oadi ng  moduli 
defined  by  the  equations 


Ks  -  Ks  (Ipk) 


G  =  G(  k) 


(R.15) 

(R.16) 
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The  particular  cap  model  discussed  here  is  that  presented  by  Baladi 
and  Akers  (1981),  except  that  the  shear  failure  criterion  has  the  equation 
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.  .  - -  ^  ^  - 


JJ2  =  f(Ij)  =  (Cs  ♦  alj)  ♦  C(1  -  e  x) 


( R .  1  / 


The  ellipsoidal,  strain  hardening  cap  yield  surface  has  the  equation 
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lh  - L ) 

2 

+ 

X  -  L 

[f  (k ) 

=  1  (L  <  Ij  £  X) 


(R.lt 


so  that 


^  =  Fdj.k)  .  f(k) 


{Ix~d2  -  (IrL)2  =  £  /r  '  '2  ‘  '2 


^  H  (X-L)  -  (Irl)‘ 


(R.l! 


where 


R 


X-L 

=  TR) 


(R.2C 


R  is  the  ratio  of  the  cap  ellipse  principal  axes,  but  it  is  not  obvious 
which  is  the  major  and  which  the  minor  principal  axis  until  the  model 
parameters  have  been  Getermined.  The  relation  between  plastic  volumetric 
strain,  e^,  and  peak  hydrostatic  effective  stress,  of 
determined  by  hydrostatic  unloading,  is 


70CT’ 


-  W 
eKK  "  w 
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e“3D(o0CT  "  GrM 


(R.21 


where  K  and  D  are  material  constants,  and  Cp  is  the  hydrostatic 
component  of  the  geostatic  effective  stress  tensor.  This  relation  is 
shown  in  Figure  (R.2).  The  relation  between  X,  the  value  of  1^  at  which 
the  ellipsoidal  cap  meets  the  Ij  axis,  is  patterned  after 
Equation  (R.21).  However,  a  mathematical  device  is  inserted  to  prevent 
further  shrinking  of  the  cap  due  to  dilatancy  on  the  shear  failure  surface 


when  the  corner  has  shrunk  to  the  axis.  This  means  k  cannot  be 
negative.  The  relation  between  X  and  is  constructed  by  defining 
an  auxiliary  plastic  volumetric  strain  variable,  z,  for  which 


dz  =  0 


( k  _<  0  and  dc^  <  0) 
(otherwi se) 


(R.22 


so  that  z  will  decrease  due  to  dilatancy  on  the  shear  failure  surface  only 


when  k  >  0.  Then  z  is  related  to  X  by  the  equation 


— 

=  /  dz  =  W  1 


e-D(X  -  3Gr) 


(R.23 


Inversion  of  Equation  (R.23)  yields 


*  ■ 3Gr  *  ®  (r1! 


(R.24 


The  parameters  W  and  D  in  Equations  (R.21),  (R.23),  and  (R.24)  can  be 
obtained  from  a  drained  isotropic  compression  test  in  which  e£K 
increases  continuously  so  that  and  z  are  identical. 

The  relations  between  X,  k  =  L,  R,  and'z  can  be  obtained  by  combining 
the  results  of  drained  isotropic  compression  tests  and  drained  triaxial 
compression  tests  at  constant  cell  pressure.  For  drained  triaxial 
compression  tests  at  constant  cell  pressure,  Equations  (K.3)  and  (K.9) 


so  that 


!1  =  3o3C  +  V3J2 


(  r  .  2  7 : 


Figure  (R.3)  shows  the  hypothetical  results  of  such  a  test. 

Figure  (R.3a)  shows  the  loading  stress  path,  with  slope  When  the 

stress  point  is  at  P,  the  flow  rule  predicts  the  plastic  strain  increment, 
deP,  will  be  normal  to  the  cap  surface  passing  through  P,  and  therefore 
de£K  >  0  at  all  points  on  the  stress  oath  until  failure  occurs  at 
Point  S  on  the  shear  failure  surface.  Only  then  does  dilatancy  suddenly 
occur,  according  to  the  theory  of  the  cap  model.  Figure  (R.3b)  shows  the 
stress-strain  response,  assuming  strain  control.  Failure  occurs  at  Point  S, 
at  which  (o^  -  o^q)  is  a  maximum.  From  unloading  data  a  curve  showing 
the  plastic  volumetric  strain,  e|^,  corresponding  to  a  given  value  of 
°1  ‘  °3c  on  the  loadln9  stress-strain  curve  can  be  constructed,  and  from 
that  curve  the  plastic  volumetric  strain  at  failure,  e£k  can  be 
found.  The  triaxa'a!  compressive  strength,  (o^  -  o^)^,  is  the 
ordinate  at  Point  S  on  the  stress-strain  curve,  so  that  Equations  (R.25)  and 


(R.26)  yield 


1  If  =  (ol  "  °3C)f  +  3a3C 


(R.28) 


(ol  "  °3C)f 


(R.29) 


Figure  (R.4)  shows  how  the  results  of  isotropic  and  triaxial 
compression  tests  can  be  combined  on  a  three-axis  plot  which  yields  X, 
k  =  L,  and  R  as  functions  of  ej^.  For  computational  purposes 
e[L  is  replaced  by  z  to  prevent  excessive  shrinking  of  the  cap.  The 
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above  discussion  is  conceptual.  Cap  model  users  often  determine  or  refine 
the  model  parameters  by  trial  and  error. 

The  drained  hypoelastic  incremental  unloading/reloading  bulk  modulus 
is  defined  by  the  equation 


K_  = 


Kis 


is 


1  - 


KlSe 


'K2S(I1 


-  3Gr) 


( R . 30 ) 


and  the  hypoelastic  incremental  unloading/reloading  shear  modulus,  which 
applies  to  both  drained  and  undrained  conditions,  is  defined  by  an 
equation  of  the  same  form 


The  parameters  Kis,  K^,  K2S,  ,  Gj,  and  G^  are  material 
constants.  The  shear  modulus  is  obtained  from  the  unloadi ng/rel oadi ng 
portion  of  a  plot  of  (oj  -  a-j)  against  (ej  -  e^),  measured  auring 
a  triaxial  compression  test  with  both  axial  and  radial  or  volumetric 
strain  measurements. 

Equations  (R.17),  (R.21),  (R.30),  and  (R.31)  are  all  of  the  form 


y 


a  +  bx 


ce 


-dx 


(R.32) 


the  parameters  of  which  can  be  determined  as  follows: 

V  1  V  1 

1.  Plot  i  against  -.  The  value  of  =£■  extrapolated  to  —  — >0  is  b. 

a  A  A  A 

See  Figure  (R.5a). 


2.  Plot  y-bx  against  i.  The  value  of  y-bx  extrapolated  to  i  >  0  is  a. 
See  Figure  ( R . 5b) . 
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3.  Plot  in  (a  +  bx  -  y)  against  x.  The  y  intercept  is  In  c  and  the 
slope  is  -d,  since  Equation  (R.32)  yields 


1 n( a  +  bx  -  y)  =  ln(ce~ax)  =  In  c  -  dx 
See  Figure  (R.5c) . 


( R .  33 


The  cap  model  response  to  prescribed  strain  inputs  is  hypoelastic, 
provided  the  stress  point  lies  within  the  yield  surfaces.  In  this  case 


Equations  (R.15)  and  { R . 1 6 )  yield 
dIl  =  3Ksd£KK 


( R .  34 


dSij  =  2G  deij 


If  the  assumed  hypoelastic  stress  increment  (the  elastic  trial) 


(R.35 


violates  the  tension  cutoff  the  first  stress  invariant  is  set  equal  to  -T, 
and  all  deviator  stresses  are  set  equal  to  zero.  In  this  case  the  plastic 
volumetric  strain  increment  is  calculated  from  the  equation 


so  that 
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(R.36 
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If  the  elastic  trial  point  lies  to  the  left  of  the  corner  and 


violates  the  shear  failure  surface,  the  stress  point  is  corrected  back  to 
the  shear  failure  surface  using  the  method  of  Appendix  P  for  a  perfectly 
plastic  material.  Since  the  cap  model  is  entirely  associative, 

Equations  (P.17),  (P.23),  and  (P.27)  yield 


The  cap  shrinks  (if  possible)  in  response  to  the  dilatant  plastic 
volunetric  strain  given  by  Equation  (R.39). 

Equations  (R.40)  and  (R.41)  define  a  point  on  the  shear  failure 
surface  (extended  beyond  the  corner,  if  necessary),  but  that  point  will 
violate  the  cap  if  the  point  lies  to  the  right  of  the  shrunken  corner. 

This  situation  is  shown  in  Figure  (R.6),  where  N  is  the  initial  stress 
point,  E  the  elastic  trial  point,  T  the  shear  failure  surface  correction 
point  assuming  the  plastic  strain  increment  to  be  directed  along  normal  PQ 
to  the  shear  failure  surface,  C  the  initial  corner,  and  D  the  shrunken 
corner  (which  may  coincide  with  C).  If  the  shear  failure  surface 
correction  point  had  been  taken  as  P  instead  of  T,  implying  a  plastic 
strain  increment  lying  along  vertical  PR,  the  corner  would  have  remained 
at  C  and  cap  violation  would  not  have  occurred.  As  the  assumed  shear 
failure  surface  correction  point  moves  from  P  toward  T  along  the  shear 
failure  surface  the  direction  of  the  plastic  strain  increment  shifts  from 
PR  toward  PQ  and  the  corner  point  shrinks  (if  it  can)  from  C  toward  D. 
Eventually  the  assumed  corrected  stress  point  meets  the  shrunken  corner 


The  derivative  dk/dej;^  is  determined  empirically.  (See 
Figure  (R.4).)  The  final  aeviator  stresses,  after  all  the  above 
corrections  have  been  applied,  are  given  by  the  equation 


(R.47) 


If  the  elastic  trial  point  lies  to  the  right  cf  the  corner  and 
violates  the  cap,  a  trial  and  error  correction  is  used  (in  place  of  the 
incremental  equations  of  Appendix  G).  The  correction  places  the  final 
stress  point  on  the  expanded  cap  and  also  satisfies  a  secant  flow  rule. 
The  expanded  yield  surface  has  the  equation 


[of1  .  Fdp.k"*1) 
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and  the  flow  rule  takes  the  form 
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Note  that 


a  =  0 


(Ix  =  L) 
(I1  =  X) 


(R.51) 


Equation  (R.49)  yields 


Eliminating  d>.  between  Equations  (R.54)  and  (R.57)  yields 


so  that 
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Equations  (R.48)  and  (R.59)  are  the  two  equations  which  the  final 
stress  point  must  satisfy. 

Figure  (R.7)  illustrates  the  trial  and  error  solution  process,  which 
proceeds  (conceptually)  as  follows: 

1.  Assume  k(=  L) . 

2.  Determine  dej^,  R,  and  X  from  empirical  relations.  [See 
Figure  (R.4)]. 

3.  Compute  1^  from  Equation  (R.54). 

4.  Compute  a  from  Equation  ( R . 50 ) . 

5.  Compute  Jj^  from  Equation  (R.59).  Steps  3  and  5  locate  a  point 
on  the  flow  path  in  Figure  (R.7). 

6.  Compute  from  Equation  (R.48).  Steps  3  and  6  locate  a  point 
on  the  yield  path  in  Figure  (R.7). 

7.  The  flow  path  and  yield  path  eventually  intersect  at 


The  flow  path  cannot  pass  to  the  left  of  the  expanded  corner  and  miss 
the  yield  surface  altogether,  because  Equations  (R.51)  and  (R.59)  show 
that  the  flow  path  would  have  a  vertical  secant  if  1^  were  equal  to  L. 
Therefore,  the  flow  path  will  always  intersect  the  yield  path  (or  the 
expanded  cap)  to  the  right  of  the  expanded  corner.  Also,  since 
Equations  (R.51)  and  (R.59)  show  that  the  flow  path  has  a  horizontal 
secant  when  1^  =  X,  it  is  convenient  to  (conceptually)  set 


a 
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so  that  the  flow  path  will  be  a  horizontal  straight  line  until  1^  <  X. 

Based  on  the  above  analysis,  it  is  apparent  that 

k°  <  k  <  l[  (R.61) 


so  that  the  range  of  k  to  be  tested  in  the  trial  and  error  solution  of 
Equations  (R.48)  and  (R.59)  is  bounded. 

Simultaneous  solution  of  Equations  (R.48)  and  (R.59)  actually 

n+  1 

consists  of  finding  the  value  of  k  =  kn  1  for  which 

kn  <  k  <  l[ 

L  <  Ij  <  X 


and 
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The  function  on  the  LHS  of  Equation  (R.63)  is  not  bounded,  so  it  is 
convenient  to  work  with  the  bounded  function 


P(k)  = 


(L  <  I,  <  X) 


( 


for  which  -1  <_  P ( k )  _<  1,  and  which  has  the  same  root  as  Equation  (R.63) 
because  the  denominator  is  always  positive.  The  function  P(k)  decreases 
steadily  as  k  increases.  8ecause  it  is  possible  to  have  Ij  >  Xn, 


it  is  desirable  to  also  define 


P(k)  = 


!1  -  L 
X  -  L 


(I:  >  X) 


(R.66) 


sc  that  P(k)  will  decrease  steadily  as  k  increases  (and  1^  decreases) 
even  if  I ^  >_  X,  and  will  be  equal  to  1  when  1^  =  X  and  therefore  never 
less  than  the  upper  bound  on  P(k)  for  the  range  (L  £  1^  _<  X).  And 
because  it  is  possible  to  have  1^  <  L  during  the  trial  and  error 


process,  it  is  desirable  to  define 


P(k)  = 


]2  -  X 


(IL  <  L) 


(R.67) 


so  that  P(k)  will  also  decrease  steadily  as  k  increases  (and  1^ 
decreases)  even  if  ^  £  L,  and  will  be  equal  to  -1  when  I ^  =  L  and 
therefore  never  greater  than  the  lower  bound  on  P(k)  for  the  range 
(L  <_  Ij  £  X).  Equations  (R.66)  and  (R.67)  ensure  that  the  final  value 
of  +  1  will  lie  in  the  interval  ln+*  <_  l"+*  £  Xn  and 

also  that  Equation  (R.60)  will  not  be  needed  because  the  quantity  a  will 
be  used  to  calculate  P(k)  only  when  1^  lies  in  the  desired  interval. 

Once  kn+1  =  Ln+1  is  found,  de£K,  Xn+*,  and  Rn+*  are 
determined  empirically,  is  calculated  by  Equation  (R.54), 


by  Equation  (R.48),  and  the  deviator  stress  components  by 
Equation  (R.47). 

When  dealing  with  finite  stress  or  strain  differences  during 
unloading/reloading,  or  during  any  phase  of  isotropic  compression,  direct 
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integration  of  Equations  (R.30)  and  (R.31)  is  sometimes  advantageous 
Both  are  of  the  general  form 


Be"Cy 


so  that 


dx  = 


eCydy 
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For  drained  hydrostatic  unloading/reloading 


Equations  (R.79)  and  (R.81)  combined  give  an  expression  tor  the 
increase  in  total  volumetric  strain  between  any  two  stress  points  on  the 
drained  loading  hydrostat.  The  inverse  calculation  can  be  accomplished 
graphically,  as  follows: 

n  0 

1.  Stipulate  p  and  ej^  ^  ^  =  0. 


2.  Assume  e 


KK,2' 


3.  Calculate  1^  2  by  Equation  (R.82). 

4.  Calculate  ej^  2  by  Equation  (R.79). 

5.  Plot  2  against  2  as  in  Figure  (R.8),  and  find 
eKK,2  as  shown. 


6.  For  a  given  e^,  enter  the  figure  backward  and  find 
6  P 

e^K ,  and  e£|p  either  of  which  can  be  used  to  find  Ip 
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Figure  R.4. 


Relations  between  cap  parameters  and  plastic 
volumetric  strain. 


(a)  Flow  path,  yield  path  and  solution  point. 
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(b)  Functions  to  be  equalized  by  trial  and  error 


Figure  R.7.  Cap  violation  correction. 


*** 


APPENDIX  S 


UNDRAINED  CAP  MODEL  AND  COMPUTATIONAL  ALGORITHM 


During  undrained  isotropic  compressive  loading  the  measured  total 
stress  incremental  bulk  modulus  is  assumed  to  have  the  same  form  as 


Equation  (R.30). 


K.  -K-  (I.-3GJ 

,/  i  rn  11/  1  p 

m  1  -  Klm  j_  lm 


(S.I) 


so  that,  by  analogy  wi th  Equations  (R.79)  and  (R.80), 
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During  undrained  isotropic  compression  Equations  (R.79),  (R.00), 
(R.81),  and  (R.82)  apply  to  the  effective  stresses,  so  that 
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The  value  of  I ' ^  ^  corresponding  to  a  given  total  volumetric  strain 
increment,  ^  ~  ekk  1*  can  calculated  using  Figure  (R.8),  and 
either  Equation  (S.5)  or  (S.7). 

For  undrained  triaxial  compression  at  constant  cell  pressure  the 
slope  of  the  total  stress  path  in  dpjd^)  space  is  1/JT,  so  that 

ll  =  3o3C  +i33 2  =  3o3C  +i^2  (S‘ 


si  nee 


The  procedure  for  calculating  an  effective  stress  path  increment  is 


as  follows: 


1.  Stipulate  ae^  and  assume  Ae^. 

I  I  I 

2.  Calculate  the  effective  stresses  o3>  Ip  and 
J J^"  using  the  drained  (effective  stress)  model. 

3.  Calculate  by  Equation  (S.3). 

4.  Calculate  the  excess  pore  pressure  by  the  equation 

..  h-1! 


5.  Calculate  the  total  radial  stress  by  the  equation 


c3  =  °3  +  u 


6.  If  a7  ±  o,r  select  a  new  value  of  ae^  and  repeat  steps  2-6. 


7.  When  o^  =  o^,  calculate  the  total  axial  stress  by  the  equation 


=  Oj  +  U 


(S.12) 


Equation  (S.8)  will  be  satisfied  because  it  takes  the  form  of 
Equation  (R.25). 

Strain-controlled  undrained  response  is  calculated  as  follows: 


1.  Stipulate  Ae... 

'  vJ 


2.  Caculate  and  Ij  using  the  drained  (effective  stress! 


model . 


3.  Calculate  1^  by  Equations  (S.3). 

4.  Calculate  the  excess  pore  pressure  by  Equation  (S.10) 

5.  Calculate  the  total  stresses  by  the  equation 


a .  .  =  o'.  +  U6  •  • 
1J  1J  U 


(S.ll) 

Notice  that  this  approach  recognizes  dilatancy  in  the  effective  stress 


model,  but  not  in  the  total  stress  model. 


APPENDIX  T 


LADE  MODEL  CROSS  SECTIONS  AND  PARAMETER  DETERMINATION 


Lade's  failure  criterion  (the  expansive  yield  surface  at  its  maximum 


extent)  has  the  equation 


f.  'I  27 1  hf 
f- =  V?i '  7  W  = 


Nov:  Equation  (C.5)  yields 


ll  =  3o0CT 


Equation  (B.24)  gives 


J1J2  ll 
h  =  J3 - 3"  +  77 


Equation  (A. 63)  yields 


/J?\3/2 

J  2  =  2  cos  3u  — 1 


and  Equation  (C.7)  yields 


Substitution  of  Equations  (T.3),  (B.24),  (T.4),  and  (T.5)  into 
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then  Equation  (T.7)  can  be  written  in  the  form 


Az3  -  1.5z2  +  B  =  0 


or 


(1.5  -  Az)  =  B 


(T.12) 
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or 
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(T.13) 


By  fixing  u  and  therefore  A  in  Equation  (T.10),  and  varying  h  and 
therefore  B  in  Equation  (T.11),  Equation  (T.13)  will  yield  a  longitudinal 
cross  section  of  the  failure  surface.  By  fixing  h  and  therefore  B,  and 
varying  u  and  therefore  A,  Equation  (T.13)  will  yield  a  transverse 
(octahedral)  cross  section.  For  triaxial  compression 
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so  that  if 
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where  (5  is  the  compression  angle  of  obliquity,  then 


ll  =  °1  +  2o3  =  °3  (N0  +  2) 
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and  therefore  Equation  (T.2)  yields 
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where  Equation  (T.ll)  yields 


B  =  1  -  27r 
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Values  of  r,  C,  and  h  corresponding  to  compression  angles  of 
obliquity  of  30,  35,  40,  45,  and  50  degrees  are  tabulated  below: 


p 

r 

B 

h 

DEG 

— 

— 

— 

30 

C. 02400 

0.35200 

0.06818 

35 

0.02003 

0.45920 

0.04362 

40 

0.01600 

0.56788 

0.02818 

45 

0.01215 

0.67199 

0.01808 

50 

0.00867 

0.76590 

0.01132 

Values  of  z  for  Lade's  failure  criterion  are  tabulated  in 
Table  (T.l).  The  rows  show  longitudinal  cross  section  values  for  a  given 
-o,  and  the  columns  show  octahedral  cross  section  values  for  a  given  h. 

The  determination  of  Lade's  model  parameters  proceeds  as  follows. 

The  unloadng/reloading  elastic  modulus,  Eur,  is  measured  on  plots  of 
^al  *  °3C^  a9dtnst:  ei  from  drained  triaxial  compression  tests  at 
constant  cell  pressure,  and  is  assumed  to  be  an  exponential  function  of 
cell  pressure  of  the  form 

EUr  '  Kur  (^f)  <E-21 

Taking  the  logarithm  of  both  sides  of  Equation  (T.20)  yields 

'°9lOEUr  =  '^loV  *  n  '°910  (Jfj  0.2 

so  that  the  parameters  Kur  and  n  can  be  obtained  from  a  log  log  plot  of 
E  against  o^f/p  ,  which  should  appear  as  a  straight  line. 


The  unloading/reloading  bulk  modulus,  Bur,  can  be  obtained  from 
drained  isotropic  compression  tests  with  unloading,  and  since 


ur  3(1  -  2vur) 


(T .22) 


where  vur  is  Poisson's  ratio  for  unloading/reloading,  the  value  of  vur 
can  be  found  from  the  expression 


V  =  7 


(T.23) 


Frequently,  vur  is  found  to  be  relatively  insensitive  to  o^,  and 
is  therefore  assumed  constant. 

Drained  isotropic  compression  activates  only  the  collapse  yield 
surface,  so  Wc  can  be  determined  from  drained  isotropic  compression 
tests  with  unloading  as  the  area  under  a  plot  of  octahedral  normal  stress 
against  plastic  volumetric  strain.  For  isotropic  compressive  primary 
loading  Equations  (3.8.1),  (3.8.2)  and  (3.8.3)  yielc 
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Taking  the  logarithm  of  both  sides  of  Equation  (T.25)  yields 
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so  that  the  parameters  C  and  p  can  be  obtained  from  a  log  log  plot  of 


W  /p,  against  3(orirT/p,)‘:,  which  should  appear  as  a  straight 

r  *  UL I  a 


c'  ra 
line. 


Shear  failure  in  a  strain-control  led  drained  triaxial  compression 
test  at  constant  cell  pressure  is  defined  as  attainment  of  the  maximum 
value  of  (oj  -  o^).  At  this  point  Equation  (3.8.8)  yields 
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so  that 
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Taking  the  logarithm  of  both  sides  of  Equation  (T.28)  yields 


,3 


lo9l0  (t~  -  2j  f  =  1O910  nl  -  m  1O910  UM  f 


(T.29) 


so  that  the  parameters  and  m  can  be  obtained  from  a  log  log  plot  of 

f1!  )\  M 

U - 27  against  —  ,  which  should  appear  as  a  straight  line. 

V3  /  f  \pa/  f 


Having  the  parameters  Eur,  vur,  C,  and  p,  the  elastic  and 


collapse  plastic  strains  can  be  computed  throughout  a  triaxial  compression 
test  at  constant  cell  pressure.  The  elastic  and  collapse  plastic  axial 
anc  volumetric  strain  increments  are 
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and  the  expansive  axial  and  volumetric  strain  increments  are 


deP  = 

dei  - 

(de*  +  dej ) 

(T 
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work  i  s 
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The  values  of  Wp 

correspond!'  ng 

1  1 
to  f  =  and  fp  =  0.60 

respectively,  are  called  Wp  pEA|<  and  Wp  60. 

The  assumption  is  that  Wp  pE^  is  related  to  o-^  by  the  expression 

VPEAK  ■  ""a  (£)'  (T 

Taking  the  logarithm  of  both  sides  of  Equation  (T.37)  yields 

.  !•.,  J  .  1  lo8ln  M  (T 


so  that  the  parameters  F  and  1  can  be  obtained  from  a  log  log  plot  of 


Up  PEAK  °3f 

against  — -,  which  should  appear  as  a  straight  line. 


Now  because  f  and  fp  are  equal  as  long  as  (o^  -  o^) 

II 

is  increasing,  and  fp  is  a  function  of  Wp  according  to 
Equation  (3.8.9),  let  Wp  P£^  be  the  value  of  Wp  which  maximizes 

I)  I 

fp  (or  fp)  at  the  value  Then  Equation  (3.8.9)  yields 
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so  that  at  failure 


and  therefore  since  is  not  zero, 


Substitution  of  Equation  (T.41)  into  Equation  (3.8.9)  yields 


so  that  according  to  the  definitions  of  Up ( peak  anci  WP,60» 


Note  that  Wp  must  be  measured  as  a  function  of  o^,  but 
empirical  formula  relating  Wp  to  is  needed.  Dividing 
Equation  (T.44)  by  Equation  ( T .43 )  yields 


0.60  = 


so  that 


P.60 

'p.PEAK 


,  WP ,60 
i  -  v — '  — 
P.PEAK 


(T  .45 ) 


'P,PEAKj 
In  0.60 


+  1  - 


P.PEAK 


(T.46) 


The  assumption  is  that  q  is  related  to  by  a  linear  equation  of 


the  form 


q  =  a  +  8 


(T.47) 


so  that  the  parameters  a  and  6  can  be  obtained  from  an  arithmetic  plot  of 
q  against  should  appear  as  a  straight  line. 

Having  both  Wp  and  q  as  functions  of  o^,  t^1e  parameter  b 

can  be  obtained  as  a  function  of  o^q  by  Equation  (T.41),  and  the 
parameter  a  can  be  obtained  as  a  function  of  by  writing 

Equation  (T.43)  in  the  form 


The  parameter  no  in  Equation  (3.8.11)  can  be  determined  from  the 
expansive  plastic  Poisson's  ratio. 


obtained  from  a  drained  triaxial  compression  test  at  constant  cell 
pressure,  where  using  the  results  of  Equations  (T.34)  and  (T.35)  yields 

dc3  =  ?(deS  '  del}  (T 

so  that 


de3  1 


(T 


Using  the  flow  rule  and  Equation  (3.8.11)  yields 


0.56989 
0.S4775 
0.89779 
0.84263 
0.79238 
C. 74967 
0.71456 
C. 68531 
0.66106 
0.64216 
0.62872 
0.62068 
0.61801 


0.87226 
0.85720 
0.82084 
0.77766 
0.73624 
0.69987 
0.66932 
C.  64356 
0.62205 
0.60520 
0.59317 
0.58596 
0.58356 


0.77128 

0.76119 

0.73557 

0.70321 

0.67057 

0.64086 

0.61530 

0.59342 

0.57501 

0.56050 

0.55007 

0.54381 

0.54173 


0.66861 

0.66199 

0.64456 

0.62138 

0.59687 

0.57374 

0.55330 

0.53551 

0.52040 

0.50840 

0.49973 

0.49450 

0.49275 


0.56569 
0.56149 
0.55012 
C. 53434 
0.51694 
0.49991 
0.48442 
0.47071 
0.45893 
0.44949 
0.44263 
0.43847 
0.43708 


APPENDIX  U 


ARA  CONIC  MODEL  CROSS  SECTIONS  AND  PARAMETER  DETERMINATION 


Triaxial  cross  sections  of  the  conic  model  compressive  and  expansive 
yield  surfaces  appear  as  shown  in  Figure  (U.l). 

The  conic  model  failure  criterion  (the  expansive  yield  surface  at  its 
maximum  extent)  has  the  equation 


fi‘ 


‘OCT1 


1  -  E  cos  3u 


Vc0CT 


+  m  =  n 


1 


(U.l) 


or 


’1 


°oct' 


l0CT 


1  -  E  cos  3o p. 

/°oct\ 

1  +  m 


(U.2) 


When  «  is  fixed  in  Equation  (U.2)  the  variation  of  tqct  with 
o0CT  becomes  hyperbolic,  taking  the  form 


l0CT 


(  nl  ^ 

|c0CT 

\  1  -  E  cos 

1  +  1 

( m  ) 

1  °0CT 

(U.3 ) 


The  initial  slope  of  the  failure  surface  at  zero  octahedral  normal 
stress  is 


(dx 


OCT 


vOCT/o  1  -  E  cos  3u 


(U.4 ) 


and  the  upper  shear  strength  asymptotic  limit  is 
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When  aQQ-p  is  fixed  in  Equation  ( U . 2 )  the  variation  of 
TQci/aQcT  becomes  a  triple  ellipse,  taking  the  form 

_ ^1 

1  +  m 

t0CT  _ 

Gqq j  1  -  E  cos  3u 


The  parameters  n^,  m,  and  E  can  be  determined  from  a  series  of 
triaxial  compression  and  extension  tests.  For  triaxial  compression 
('-  =  120°),  Equation  (U.3)  reduces  to 


or 


aOCT  _  1  -  E 
r0CT  nl 


where 


k 


lc 


1  -  E 


nl 


For  triaxial  extension  (<j  =  60° ),  Equation  (U.3)  reduces  to 


Plots  of  Equations  (U.7)  and  (U.10)  are  shown  in  Figure  (U.2). 


Having  determined  the  parameters  k^c,  k^c  >  klp,  and  k^,  the 
parameter  m  can  be  calculated  from  the  expression 


m  = 


2c 


2e 


rr  =  ft 


lc 


le 


which  provides  a  consistency  check.  The  parameters  and  E  can  be 
calculated  from  Equations  (U.8)  and  {U.ll),  written  in  the  form 

kic-i  +  E  -  1 

k le° 1  -  E  =  1 


so  that 


2k 


klcril  =  1 


1  c 


1(  +  1(  * 
Kle  Klc 


kle  '  klc 
kle  +  klc 


(II.  15) 


Having  calculated  the  parameters  and  m  (and  E)  from  triaxial 
compression  and  extension  tests,  the  accuracy  of  the  assumed  octahedral 
cross  section  form  can  be  investigated  by  a  series  of  true  triaxial 
tests.  If,  in  Equation  (U.6)  we  set 


’1 


1  +  m 


ff0C  T 


(U. 16) 


then  Equation  (U.6)  can  be  written  in  the  form 

T0CT  _  4> 

Oq^  1  “  E  cos 


(U .  17 ) 


which  is  the  equation  of  a  triple  ellipse  in  polar  coordinates. 

Equation  (U.17)  can  be  written  in  a  linear  form  to  obtain  the  octahedral 
eccentricity,  E,  as  a  consistency  check  on  the  previously  determined  value 
from  Equation  ( U.15) . 


T 


OCT 

OCT 


1  -  E  cos  3« 


(U. 18) 


A  plot  of  Equation  (U.18)  is  shown  in  Figure  (U.3). 

The  method  for  computing  the  conic  model  octahedral  cross  section  sets 
U£  =  <J  (A. 65b) 

in  Figure  (C.2),  so  that 

tan  u  =  - -  ( U.  19) 

u 
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where  u  is  Lode's  parameter.  And  if  we  assume  that  d  =  0  in  Figure  (C.2) 
and  Equation  (C.26),  then 


O'P'  =  ^3  -2£I  =  z  (U.2 

°0CT 


so  that  the  horizontal  projection  of  O'P'  can  be  calculated  in  two  ways: 


‘3*  z  sin  u  =  ^6  -  z  cos  uj  (U.2 


Solving  Equation  (U.21)  for  z  yields 

z  (^3  sin  <0  +  cos  u  sin  =  /?  sin  0 


or 


z  = 


v2  sin  0 


S3  sin  u  +  cos  u  sin 


and  solving  Equation  (U.21)  for  sin  0  yields 
y/3  z  sir,  u 


sin  0  = 


7  -  Z  COS  u 


Note  that  when  m  =  0,  Equations  (U.7)  and  (U.10)  yield 


1 


kic  =z^ 


1 


k.  =  — 
lfc  ze 


(U.2 


(U.2 


(U.2 


(U.2 


Octahedral  cross  section  data  for  the  case  (m  =  0;  0C  =  32°; 

=  35°)  are  tabulated  in  Table  U.l  and  plotted  in  Figure  (U.4).  The 
calculation  sequence  used  to  obtain  the  values  shown  in  Table  U.l  and 
Figure  (U.4)  is  as  follows: 


M 


(U.19) 


2  ft  sin 
c  3  -  s  i  n  # 

2  ft  sin 
e  3  +  s  i  n 


2z  z 
c  e 

z  +  z 


(«  =  120°)  ( U.26) 


(u  =  60°)  (U.27) 


(U.28) 


E 


(U.29) 


z  = 


"1 _ 

E  cos  3c 


sin  . 


ft  z  sinu 
2  -  Z  COSu 


(U.30 ) 

(U.31 ) 


De term' nation  of  the  unloading/reloading  Young's  modulus  parameters, 
«ur  and  n  in  Equation  (4.6),  for  the  conic  model  is  accomplished  in  the 
same  way  as  for  the  Lade  model,  as  described  in  Appendix  T.  The  same  is 
true  for  the  unloading/reloading  Poisson's  ratio,  vur. 

At  the  present  time  the  parameters  A,  M,  \,  y,  and  8  in 
Equation  (4.10)  are  determined  by  trial  and  error.  When  the 
unloaaing/reloading  hysteresis  option  has  been  fully  implemented,  a  method 
for  determining  the  parameters  will  be  developed. 

The  parameters  C  and  p  in  the  compressive  hardening  rule, 


Equation  (3.8.3),  can  be  determined  in  the  same  way  as  for  the  Lade  model, 
as  described  in  Appendix  T,  but  there  is  an  easier  way  which  avoids  the 


numerical  integration  required  to  obtain  Wc  as  a  function  of  fc.  It 


involves  directly  fitting  the  plastic  hydrostat  (the  curve  of  oqct 


against  obtained  from  drained  isotropic  compression  tests  with 


unloading).  Equation  (T.25)  can  be  written  in  the  form 


Wc  =  Cpa 


^°0CT ' 


( U. 32 ) 


so  that  the  compressive  plastic  work  increment  generated  by  a  compressive 
plastic  volumetric  strain  increment  is 


dWc  =  c0CTdev  =  Cppa 


'°0CT 


P-1 


.  6 


OCT 


do 


OCT 


=  S6Cp 


p-1  do 


OCT 


OCT 


and  therefore 


tfcj  =  6Cp 


k2 


p-1  /do 


0CT\ 


Integration  of  Equation  (U.34)  yields 


.c  /  6 Cp_\  Lp-1)  f°0CT)  2p ,  f°0Cj)  0 


v  =  \2p  -  1 


=  L 


where 


L .  c  l  A) 


Q  =  2p  -  1 

Taking  the  logarithm  of  both  sides  of  Equation  (U.35)  yields 


(U.33) 


(U.34) 


(U.35) 


(U.36) 

(U.37) 


n 


C. 
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1  og 


10 


1og10  L  ♦  0  log1Q 


(U. 38) 


so  that  the  parameters  L  and  Q  can  be  obtained  from  a  log  log  plot  of 

Q 

ev  against  ^OCT/Pa*  which  should  appear  as  a  straight  line.  The 
parameter  p  can  then  be  determined  by  writing  Equation  ( 1) . 3 7 )  in  the  form 


P 


Q  +  1 
2 


(U.39) 


and  the  parameter  C  can  then  be  determined  by  writing  Equation  (11.36)  in 
the  form 

c  =  (2Vi)(3l"p)  L  ,u',0) 

The  parameter  r  in  Equation  (4.3)  is  not  determined  uniquely,  but 
rather  adjusted  by  trial  and  error  so  the  sum  of  the  calculated  elastic  and 
compressive  plastic  volumetric  strains  in  a  triaxial  compression  test  at 
constant  cell  pressure  always  exceeds  the  measured  total  volumetric  strain, 
and  the  difference  steadily  increases.  The  difference  is  minus  the 
expansive  plastic  volumetric  strain,  and  if  the  above  relationship  is  not 
maintained  it  means  the  expansive  plastic  potential  surface  is  generating 
compressive  plastic  volumetric  strain,  which  is  impossible.  The 
relationship  to  be  maintained  is 

d£*  ♦  dej  >  dEy  (U.41) 


or 


( U.42) 
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/  l-2v 


6a 


Determination  of  the  expansive  yield  surface  parameters  E,  m,  and 
in  Equations  (4.4)  and  (3.8.8)  has  already  been  discussed. 

The  parameters  a,  b,  and  q  in  the  expansive  hardening  rule, 

Equation  (3.8.9).  can  be  determined  in  the  same  way  as  for  the  Lade  model, 
as  described  in  Appendix  T,  but  there  is  anottier  way  which  uses  the  entire 

l 

curve  of  f  against  W  ,  rather  than  just  the  two  values  Wp 
and  Wp  gQ.  However,  the  method  does  involve  a  derivative. 

Equation  (3.8.9)  gives 


C  ,  P  P 

f  =  ae  K  — - 


(3.8. 


so  that 


/w  \  ,  w  ’ 

In  f '  =  In  a  -  (bp  )  —  +  -  In  — 

rx  r  a  1  n  /  n  n 


( U.  43 


Di f ferentiation  of  Equation  ( U . 4 3 )  with  respect  to  Wp/pa  yields 


f1  ,/W  ' 
P  d  P 


=  -bp 


1  f^d) 


a  q  W 


a  u  i 

iT  =  qbpa  +  q  T" 


dU 


(U.44 


The  parameters  q  and  b  can  be  obtained  from  an  arithmetic  plot  of 


p,/Wn  against 
a  p 


d 


which  should  appear  as  a  straight  line. 


The  intercept  is 


qbpa 


Pa 

WP,PEAK 


(U.4! 


and  the  slope  is  q.  Having  b  and  q,  the  parameter  a  can  be  obtained  from 

-  b  W  ,, 

an  arithmetic  plot  cf  f  against  e  p  ^p/Pa'  >  which 

should  appear  as  a  straight  line  through  the  origin  with  slope  a  as 

predicted  by  Equation  (3.8.9).  The  linear  relation  between  o,r/p,  and 

a 

q  indicated  by  Equation  (T.47),  plus  Equations  ( T . 4 1 )  and  (T.48)  apply  to 
the  conic  model  as  thc-y  oo  to  the  Lade  model. 

The  parameter  r0  in  Equation  (4.5)  can  be  obtained  from  the 
expansive  plastic  Poisson's  ratio 


V 


p 


(T.4! 


obtained  from  a  drained  triaxial  compression  test  at  constant  cell 
pressure,  as  it  is  for  the  Lade  model.  The  resulting  expression  is 


n2 


P  (vP  -  0.5) 
H  (1  +  vp) 


(T.5< 


where 


P 


E) 


(T.5: 


H  = 


1  ♦  r.i 


a0CTp  * 


The  linear  relation  between Jc jq/P# »  fp,  and 
Equation  (T.58)  applies  to  the  conic  moael  as  it 


(T 


5; 


r\2  indicated  by 
does  to  the  Lade  model. 


TABLE  U.l 


ARA  MODEL  FAILURE  SURFACE  OCTAHEDRAL  CROSS  SECTION 
FOR  (m  =  0;  =  32  DEGREES;  0e  =  35  DEGREES) 


l a 

DEG 

z 

sin  0 

-1.0 

120.000 

0.60679 

0.530 

-0.9 

117.457 

0.60589 

0.550 

-0.8 

114.791 

0.60304 

0.569 

-0.7 

112.006 

0.59810 

0.586 

-0.6 

109.107 

0.59105 

0.602 

-0.5 

106.102 

0.58198 

0.614 

-0.4 

103.004 

C. 57118 

0.625 

-0.3 

99.826 

0.55901 

0.632 

-0.2 

96.587 

0.54597 

0.636 

-0.1 

93.304 

0.53259 

0.637 

0 

90.000 

0.51938 

0.636 

0.1 

86.696 

0.50681 

0.633 

0.2 

83.413 

0.49526 

0.628 

0.3 

80.174 

0.48500 

0.622 

0.4 

76.996 

0.47620 

0.615 

0.5 

73.898 

0.46894 

0.608 

0.6 

70.893 

0.46321 

0.60C 

0.7 

67.994 

0. 45897 

0.593 

0.8 

65.209 

0.45610 

0.586 

G.9 

62.543 

0.45449 

0.580 

1.0 

60.000 

0.45398 

0.574 

Figure  U.l.  Conic  model  yield  surfaces. 


Figure  U.4.  Conic  model  failure  surface  octahedral  cross-section  for 
(m  =  0,  <pQ  =  32  degrees;  <j>e  =  35  degrees). 
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APPENDIX  V 


EVALUATION  OF  EXISTING  MODELS 

V. 1  Li  near  El  asti  c 

V.I.l  Motivation 

The  principal  advantages  of  a  linear  elastic  constitutive  model  for 
transient  soil  dynamics  problems  are  the  availability  of  or  possibility  of 
obtaining  closed  form  analytic  solutions,  and  the  existence  of  proven, 
stable  numerical  solution  techniques  for  use  when  a  closed  form  solution 
does  not  exist  and  is  too  difficult  to  obtain.  The  principal  disadvantage 
is  that  real  soil  departs  from  linear  elastic  behavior  even  at  strains  of 
the  order  of  10 

V. 1 .2  Assumptions 

The  classic  theory  of  linear  elasticity  assumes  a  material  to  be 
homogeneous  (having  the  same  stress-strain  characteri sties  at  each  point), 
and  isotropic  (having  its  stress-strain  characteristics  the  same  in  all 
directions).  The  term  elastic  implies  complete  strain  recovery  upon 
unloading,  no  matter  how  large  or  small  the  applied  stresses.  The 
stresses  at  a  point  are  assumed  to  depend  only  on  the  strains  at  that 
point,  through  a  system  of  homogeneous  linear  equations  [Timoshenko  and 
Goodier  (1970:1)].  Thus  the  linear  elastic  model  is  rate  independent. 
V.1.3  Basic  Equations 

The  homogeneous  linear  equations  relating  strains  to  stresses  in  the 
theory  of  linear  elasticity  are  of  the  form 


0 


0 


0 


o  Ko  1 


0  0  1 


0  0 


0  0 


where  the  column  vectors,  |a j  and  jej,  are  defined  by  Equations  (E.5)  and 
(E.24),  and  the  elastic  constants,  M  (the  constrained  modulus)  and  Kq 
(the  coefficient  of  lateral  stress  at  rest!,  are  defined  by  Equations 
(J.33)  and  (J.34). 

V.1.4  Parameter  Determination 

There  are  a  number  of  methods  for  measuring  the  elastic  parameters,  M 
and  K  The  constrained  modulus,  M,  can  be  measured  direct?y  in  the 
laboratory  by  a  constrained  compression  (oedometer)  test,  in  which  all 
strains  in  Equation  (V.1.1)  are  zero  except  ej.  In  a  laboratory 
hydrostatic  compression  test 


E,  -  E  o  -  Ei  = 


(V.1.2) 


al  =  °2  =  a3  =  a0CT 


(V  .1.3) 


so  that  Equation  (v.1.1)  yields 


E 


where  the  elastic  bulk  modulus,  B,  is  given  by  the  expression 


/I  +  2K  \ 

B  =  M  - 3— -  (V.1.5) 

and  therefore 

K0-  (V.1.6) 

Laboratory  test  specimens  often  suffer  from  disturbance  due  to  field 
sampling  and  specimen  preparation,  and  may  not  be  representati ve  of  an 
entire  soil  mass.  A  popular  method  for  determining  elastic  parameters 
which  are  representati ve  of  an  entire  soil  mass  is  by  field  measurements 
of  the  two  elastic  body  wave  velocities  TTerzaghi  (1043:453)].  They  are: 


'  /IT. 

'1  y  p 


compressional  wavespeed 


and 


C2  = 


IG  fM(1  -  Ko> 


=  shear  wavespeed 


where  G  =  M 


(1  -  K 


o\  =  elastic  shear  modulus 


mass  density 


(V  .1.7) 


(V .1.8) 

(V  .1.9) 


Then  the  elastic  constants,  M  and  K0,  are  given  by  the  expressions 


M 


(V.1.10) 


and 


(V.  1.11) 


V.1.5  Computed  Behavior 

The  parameters  chosen  to  represent  CARES-DRY  Sand  are  listed  in 
Table  V.1.1.  Since  the  linear  elastic  model  cannot  represent  the  behavior 
of  a  highly  non-linear,  non-elastic  material  such  as  this  over  any 
significant  range  of  stress  or  strain,  the  choice  of  these  parameters  was 
completely  arbitrary.  Computed  behavior  for  this  model  is  very 
straightforward,  as  examination  of  Figures  (V.1.1)  through  (V.1.41)  shows. 
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V. 2  Linear  Viscoelastic 


V.2.1  Motivation 


The  principal  advantages  of  a  linear  viscoelastic  constitutive  model 


for  transient  soil  dynamics  problems  are  the  availability  of  or 


possibility  of  obtaining  closed  form  solutions,  the  existence  of  proven, 


stable  numerical  solution  techniques  for  use  when  a  closed  form  solution 


does  not  exist  and  is  too  difficult  to  ohtain,  and  the  ability  of  the 


model  to  dissipate  energy.  The  principal  disadvantage  is  that  the  stress- 


strain  behavior  of  real  soil  often  departs  significantly  from  linear 


viscoelasticity. 


V.2.2  Assumptions 


The  classic  theory  of  linear  viscoelasticity  assumes  a  material  to  be 


homogeneous  and  isotropic.  The  stresses  and  stress  time  derivatives  at  a 


point  are  assumed  to  depend  only  on  the  strains  and  strain  time 


derivatives  at  that  point,  through  a  system  of  linear  differential 


equations. 


V.2.3  Basic  Equations 


The  system  of  linear  differential  equations  relating  stress,  strain 


and  their  time  derivatives  at  a  point  yields  a  set  of  hereditary  stress- 


strain  relations  of  the  same  general  form  as  Equation  (V.1.1). 


1  K  K  0  0  0 

o  o 


1  Kn 
0  0 


0  0 


K  K  1  0  0 

o  o 


(V.2.1) 


000  1-K  0 

o 


0  0  0  0  1-K  0 


0  0  0  0  0  1-K. 
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except  that  the  quantities  M  and  Kq  are  influence  functions  of  time, 
rather  than  constants,  and  the  notation  x*dy  indicates  a  superposition 


inteoral  fFung  (1%F:414)]  of  the  form 


My  =  f  x(t-T)  4^  dT 


V.2.4  Parameter  Determination 


(V.?.?) 


There  are  a  number  of  methods  for  measuring  the  parameters  which 
define  the  viscoelastic  influence  functions  M(t)  and  Kq( t )  in  Equation 
( V.?.l).  The  constrained  modulus  influence  function,  M(t),  can  be 
measured  directly  in  the  laboratory  by  a  constrained  compression 
(oedometer)  relaxation  test,  in  which  all  strains  in  Equation  (V.2.1)  are 
zero  except  tj,  and  is  a  step  function.  Similarly,  the  bulk 
modulus  influence  function. 


B(t)  =  M( t ] 


1  ♦  2K  (t) 


f V. 2.3 ) 


I  4 

■  v.w! 

•  v  •*_v. 


can  be  measured  directly  in  the  laboratory  by  a  hydrostatic  compression 


relaxation  test,  in  which  the  volumetric  strain  is  a  step  function.  It  is 
often  assumed  that  the  bulk  modulus  influence  function,  B(t),  and  the 
shear  modulus  influence  function 


G(t)  .  3[M( t)  -  3( t) ] 


(V.2.4) 


**•  • 


can  be  represented  by  simple  viscoelastic  models,  such  as  the  three 
element  model  discussed  in  Appendix  L.  The  model  parameters  defining  B(t) 
and  G(t)  can  be  determined  from  Figure  (L.2)  and  Equations  (L.30),  (L.31) 


o 


and  (L.3 2).  The  influence  function  M(t)Ko(t)  can  then  be  assembled 
using  the  expression 


M(t)K  (t)  = 


(V.2.5) 
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Table  (V.2.1)  lists  the  parameters  used  for  demonstrating  the 
three-element  viscoelastic  model  discussed  in  Appendix  L.  The  parameters 
were  chosen  to  show  a  range  cf  behavior  when  the  applied  constant  strain 
rate  varied  between  100/sec  and  0.1/sec.  Both  static  and  dynamic  uniaxial 
strain  tests  were  actually  performed  on  CARES-DRY  sand  and  reported  in 
Cargile  (1984).  The  static  tests  were  run  at  strain  rates  on  the  order  of 
0.001/sec.  The  dynamic  tests  were  run  at  strain  rates  in  the  range  of 
20-85/sec.  Some  loading-rate  effect  (stiffening)  was  observed  below 
80  MPa  [Cargile  (1984:37)]. 

V.2.5  Computed  Behavior 

Only  the  uniaxial  strain  compression  test  will  be  shown  for  this 

model.  Since  the  response  would  be  linear  and  elastic  but  for  the 

time-dependent  and  energy  damping  dashpot,  this  test  will  suffice  to 

illustrate  the  effect  of  loading  rate.  Uniaxial  strain  load-unload  cycles 

were  applied  at  four  constant  strain  rates:  100,  10,  1  and  0.1  per 

second.  Figure  (V.2.1)  shows  the  axial  stress-strain  response  for  each 

strain  rate.  Note  the  stiffening  effect  with  increasing  rate  of  load 

application  and  the  non-linear  response,  both  a  result  of  the  dashpot 

elements.  Response  of  the  bulk  element  is  shown  in  Figure  (V.2.2)  and  the 

shear  element  response  is  shown  in  Figure  (V.2.3).  At  very  high  strain 

rates,  the  element  stiffnesses  collapse  to  the  stand-alone  spring 

stiffnesses,  K^.  At  near-static  strain  rates,  the  net  stiffness  is  that 

of  two  springs  in  series: 

u  K  K 

^long-term  =  1  £ 


K,  +  K„ 


Ml  O  C  \ 


At  intermediate  strain  rates,  there  is  a  gradual  transition  between  the 


TABLE  V.2.1.  VISCOELASTIC  MODEL  PARAMETERS  FOR  CARES-DRY  SAND 


Parameter 

Symbol 

Variable 

Value 

Uni  ts 

Bulk  Spring  1 
Sti ffness 

*1 

BULK  1 

480 

x  106 

Pa 

Bulk  Spring  2 
Sti ffness 

k2 

BULK2 

300 

LO 

O 

r— 4 

X 

Pa 

Bulk  Dashpot 
Constant 

n2 

C2 

10 

x  106 

Pa-s 

Shear  Spri ng  1 
Sti ffness 

SHEAR1 

360 

x  106 

Pa 

Shear  Spri ng  2 
Sti ffness 

c,2 

SHEAR2 

114 

X 

O 

CTY 

Pa 

ns2 


Shear  Dashpot 
Constant 


CS2 


10  x  106 
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V.3  Hyperbolic 

V.3.1  Motivation 

The  hyperbolic  model  consists  of  formulas  for  computing  the  tangent 
Young's  modulus  and  Poisson's  ratio,  for  use  in  an  incremental  elastic 
analysis.  It  is  a  simple,  practical  procedure  for  representing  the 
nonlinear,  stress-dependent,  inelastic  stress-strain  behavior  of  soils 
[Duncan  and  Chang  (1970:1650)].  Values  of  required  parameters  can  be 
derived  from  the  results  of  standard  laboratory  triaxial  tests,  or  from 
more  sophisticated  test  results  if  available.  The  hyperbolic  model's 
principal  drawbacks  are  that  it  does  not  fully  account  for  stress  path 
effects  on  strength,  stiffness,  or  dilatancy. 

V.3.2  Assumptions 

The  material  is  assumed  to  be  homogeneous  and  isotropic,  and  the 
relation  between  the  major  and  minor  principal  effective  stresses  and  the 
tangent  Young's  modulus  measured  in  a  triaxial  compression  test  at 
constant  cell  pressure  is  assumed  to  hold  for  any  stress  path.  In  their 
original  hyperbolic  model  Duncan  and  Chang  assumed  a  constant  Poisson's 
ratio.  An  alternate  version  obtains  a  tangent  Poisson's  ratio  by  assuming 
a  hyperbolic  relation  between  the  major  and  minor  principal  strains, 
together  with  the  same  hyperbolic  relation  between  major  principal  strain 
and  principal  stress  difference  used  to  obtain  the  tangent  Young's 
modulus.  The  tangent  Poisson's  ratio  formulation  is  intended  only  for 
stress  levels  well  below  failure. 


*  '  **  V  -4  m  w  *  **  A**  »  "  »  ' 
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V. 3.3  Basic  Equations 


The  equation  for  the  tangent  Young's  modulus  is 


Er  =  Kpa  he 


a  Vp 


1  - 


R^(l  -  sin  fa )(a^  -  ) 


.1 2 


2(c  cos  fa  +  03  sin  fa) 


(V.3.1) 


where  =  atmospheric  pressure 


K,  n,  Rf  =  dimensionless  empirical  constants 
c,  fa  =  Mohr-Coulomb  strength  parameters 

cj,  =  major  and  minor  principal  effective  stresses 

Equation  (V.3.1)  is  derived  in  Appendix  M,  and  is  identical  to  Equation 

(M.21). 

The  equation  for  the  tangent  Poisson's  ratio  in  the  alternate  version 
i  s 


G  -  F  log 


10  pi 


VT  = 


(1  -  A) 


T 


( V.3.2) 


where 


A  = 


d(oi  -  o3) 

/°  3Y1 

Rf(l  -  sin  ) ( CTj  -  ) 

v  a/ 

2(c  COS  0  +  a 3  Sl‘n  0) 

( V.3.3) 


and,  in  addition  to  the  quantities  appearing  in  Equation  (V.3.1), 
d,  G,  F  =  dimensionless  empirical  constants 
Equations  (V.3.2)  and  (V.3.3)  are  derived  in  Appendix  N,  and  are  identical 
to  Equations  (N.7)  and  (N.8),  respectively. 


V.3.4  Parameter  Determination 

The  determination  of  the  parameters  entering  Equations  (V.3.1), 
(V.3.2),  and  (V.3.3)  is  explained  in  Appendices  M  and  N,  where  the 
equations  are  derived. 

Parameters  for  a  modified  hyperbolic  model  were  determined  for 
CARES-DRY  Sand.  The  principal  differences  between  the  model  described  in 
appendices  M  and  N  and  the  model  actually  used  are: 

(a)  Maximum  past  axial  strain  was  used  as  a  state  variable  to  define 


the  unloading  condition.  Upon  unloading,  i.e.. 


ea  <  ea,  max 


(V.3.4) 


the  tangent  Young's  Modulus  and  Poisson's  ratio  are  set  to 


constant  values, 


ET  =  Eur  =  Kur  *  p£ 


(V.3.5) 

(V.3.6) 


vT  =  vur 

(b)  Shear  stresses  were  further  constrained  to  lie  within  the 
failure  envelope  defined  by  the  Mohr-Coulomb  strength 
parameters,  c  and  t>.  Plastic  strains  were  computed  using  a 
non-associated,  Von  Mises  plastic  potential. 

(c)  Pressure  was  not  allowed  to  achieve  negative  values.  If  tensile 
failure  occured,  all  stresses  were  set  to  zero. 

Parameters  for  unloading-reloading  behavior  (Kur,  vur)  and 
strength  parameters  (c,  t>)  were  determined  by  hand  from  examination  of  the 
data  for  remolded  specimens  as  presented  in  [Cargile  (1984)].  Ail  other 
input  parameters  were  automatically  determined  in  the  SEM.  Figure 
(V.3.1)  shows  the  triaxial  stress-strain  data  used  for  this  fitting 
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process.  In  Figure  (V.3.2)  the  axial  data  has  been  replotted  on  a  Kondner 
plot.  Note  that  the  data  curves  are  very  nearly  hyperbolic,  with  most 
deviation  occurring  at  low  strains.  Figure  (V.3.3)  shows  a  logarithmic 
plot  of  Fi  vs.  a3c,  which  is  used  for  determining  K  and  n.  There  is 
considerable  scatter  in  this  fit.  Tangent  values  of  Poisson's  ratio  are 
plotted  in  Figure  (V.3.4)  against  radial  strain.  The  data  shows  wide 
variation,  and  Poisson's  ratio  exceeds  0.5  when  the  specimens  are 
dilating.  A  logarithmic  fit  for  variation  of  initial  Poisson's  ratio  with 
confining  pressure  is  shown  in  Figure  (V.3.5).  In  the  actual  model 
computation,  the  tangent  Poisson’s  ratio  is  not  allowed  to  exceed  0.5. 

V.3.5  Calculated  Behavior 

Isotropic  compression  behavior  for  the  hyperbolic  model  is  compared 
with  test  data  in  Figure  (V.3.6).  Since  triaxial  data  and  not  IC  data 
were  used  to  fit  the  loading  modulus,  and  because  the  loading  modulus  is 
strictly  dependent  on  initial  confining  pressure  and  axial  strain, 
calculated  and  actual  behavior  are  not  expected  to  match  for  this  test. 

CTC  behavior  is  covered  in  Figures  (V.3.7)  through  (V.3.10), 

Stress-strain  data,  used  for  parameter  fitting,  is  matched  well  [Figure 
(V.3.7)].  The  inclusion  of  special  unloading- reloading  modifications 
allows  the  model  to  undergo  permanent  compaction  [Figures  (V.3.7), 

(V.3.9),  and  (V.3.10)].  Since  Poisson's  ratio  is  artificially  held  to  be 
less  than  0.5,  dilation  is  not  predicted  [Figure  (V.3.9)].  Initial  bulk 
stiffnesses  are  well  matched,  as  seen  in  Figure  (V.3.10),  but  once 
dilation  commences,  calculated  and  actual  pressure-vol ume  responses 
diverge. 

Since  the  CTE  test  involves  axial  expansion,  the  model  responds  by 
unloading  with  a  high  modulus  and  the  unload-reload  value  of  Poisson's 


ratio  [Figures  (V.3.11)  and  (V.3.12)].  Presence  of  the  Mohr-Coulonb  shear 
failure  envelope  limits  stress  difference  but  does  not  predict  the  lower 
observed  values  of  shear  strength  in  extension  compared  with  those  in 
compression.  Figure  (V.3.13)  shows  that  the  predicted  volume  strains  are 
expansive  while  those  measured  were  compressive.  Pressure-volume  behavior 
is  shown  in  Figure  (V.3.14). 

The  PT C/E  and  PSC/E  [Figures  (V.3. 15-V.3. 18)  and  f V. 3. 19-V.3.20) , 
respectively]  highlight  features  of  this  model,  as  it  is  currently 
impl emented: 

(a)  Loading  and  unloading  behaviors,  as  determined  by  axial  strain, 
are  substantially  different  [Figure  (V.3.15)  for  RTC/E  and 
Figure  (V.3.19)  for  PSC/E].  This  is  due  to  the  modified  nature 
of  the  unloading  part  of  the  model. 

(b)  Without  recognition  of  the  Mohr-Coulomb  failure  surface  in  the 
model  itself,  all  loading  stress  paths  would  have  an  ultimate 
shear  strength  equal  to  (o^  -  o.jc)yLT  as  predicted  by  the 
hyperbolic  model  for  CTC  behavior.  This  would  severely 
overestimate  strength  for  certain  stress  paths  such  as  CTE  or 
RTC.  When  the  failure  envelope  is  used,  initial  stress-strai n 

-  resoonse  is  identical  to  CTC,  but  the  material  fails  in  shear 
prior  to  (oj  .  ojc)yyj.  This  can  be  seen  in  Figures 
(V.3.15)  and  (V.3.16)  for  RTC  and  in  Figures  (V.3.19)  and 
(V.3.20)  for  PSC  and  PSE. 

(c)  No  dilation  behavior  is  predicted  in  this  model  [Figure 
(V.3.17)],  and  the  pressure-vol ume  response  is  a  function  of 
initial  confining  pressure  and  axial  strain  [Figure  (V.3.18)]. 
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The  behavior  of  the  modified  Duncan  and  Chang  hyperbolic  model  is 
summarized  in  Figures  (V.3.21)  and  (V.3.22)  for  stress  paths  run  in  the 
triaxial  cell  from  7.1  MPa  initial  confining  pressure. 

Uniaxial  compression  response,  starting  at  a  low  confining  pressure 
and  tested  to  a  relatively  high  level  of  stress,  cannot  be  well  predicted 
by  this  model  as  is  shown  in  Figures  (V.3.23)  through  (V.3.25).  Axial 
stress-strain  response  [Figure  (V.3.23)]  is  much  too  stiff  because  of  the 
large  axial  strains  and  modulus  dependence  on  axial  strain.  Poisson's 
ratio  also  is  affected  and  quickly  reaches  the  0.5  limit  imposed.  This 
limit  affects  the  shear  modulus,  as  can  be  seen  in  Figure  (V.3.24). 
Unload-reload  behavior  and  compaction  are  reasonably  well  modeled. 

Predicted  uniaxial  extension  behavior  [Figures  (V.3.26)  through 
(V.3.28)]  is  essentially  the  same  as  that  described  for  the 
elastic-plastic  model  in  Section  V.5.  This  is  a  result  of  the  linear 
elastic  unload-reload  response  with  a  Mohr-Coulomb  failure  surface. 

The  axisymmetric  strain  path  exercises  disclose  some  of  the 
consequences  when  this  type  of  model  is  used  for  more  general  conditions 
like  those  occurring  in-si tu.  Predicted  behavior  for  the  WES  strain  paths 
is  compared  with  the  data  in  Figure  (V.3.29).  For  this  strain  path,  axial 
strain  is  increasing  until  the  very  last,  while  volume  strain  reverses. 
Therefore  this  model,  which  uses  axial  strain  to  signal  unloading,  never 
unloads  and  the  result  is  the  pressure-volume  response  shown.  The  bulk 
modulus  does  not  change  upon  increased  radial  expansion,  but  continues  to 
gradually  decrease  with  increasing  axial  strain.  Shear  response  is 
reasonable  but  does  not  respond  to  the  strain  path  changes  as  does  the 


The  Lade  strain  path  1  results  [Figure  (V.3.30)]  show  an  initial 
break  in  the  strain  path  and  bulk  modulus  when  the  limit  of  0.5  for 
Poisson's  ratio  is  reached.  This  results  in  a  large  departure  from 
observed  behavior.  Strain  path  2  [Figure  (V.3.31)]  also  causes  the  model 
to  stiffen  just  prior  to  the  reversal  in  volume  strain.  This  causes  a 
sharp  drop  in  pressure  (without  unloading  as  defined  by  axial  strain)  and 
the  stress  path  encounters  the  failure  surface.  Again,  deficiency  in  the 
bulk  response  causes  bad  comparison  with  the  data. 

Results  for  the  true  triaxial  strain  path  are  shown  in  Figures 
( V . 3.32 )  and  (V.3.33).  The  axisymmetric  portion  of  the  response  looks 
good,  but  when  the  strain  path  causes  stress  response  out  of  the  triaxial 
plane,  the  results  are  not  good.  Note  that  this  model  keys  on  axial  and 
radial  strains/stresses  while  ignoring  the  intermediate  principal 
strain/stresses.  Figure  (V.3.33)  shows  that  although  provision  is  made 
for  tensile  failure,  the  material  reloads  before  re-compacting. 
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TABLE  V. 3. 1(a).  HYPERBOLIC  MODEL  PARAMETERS  FOR  CARES-DRY  SAND 


Parameter  Symbol 

Variable 

Value  Units 

K 

HK 

*79.1 

Constants 

*ur 

HKUR 

4?320 

-- 

n 

HN 

0.6057 

Rf 

HRF 

0.8475 

Mohr-Coul omb 

c 

HC 

5.0  x  105 

Pa 

Stre ng th 

t 

HPHI 

30.0 

degrees 

Poisson's  Ratio 

G 

HG 

0.4924 

-- 

Constants 

F 

HF 

-0.0935 

Unload-Reload 

d 

HD 

2.004 

Poisson' s  Ratio 

Uur 

HNUR 

0.200 

Confining  Pressure 
Initial  Young's 

a3c 

HSIGMA3 

See  (b)  below 

Pa 

Modul us 

Stress  Difference 

Ei 

HE  I 

See  (b)  below 

Pa 

at  CTC  Failure 
Initial  Poisson's 

(°1- 

a3c)f  HFSDIFF 

See  (b)  below 

Pa 

Ratio 

vi 

HNUI 

See  (b)  below 

-  - 

Mass  Density 

p 

RHOREF 

1900 

kg/m3 

TABLE  V. 3. 1(b)  HYPERBOLIC  MODEL  CONFINING  PRESSURE  DEPENDENT  PARAMETERS 


a3c  E-j  HFSDIFF  v,- 

(MPa)  (MPa)  (MPa)  (  — ) 


1AF0SR  SOIL  ELEMENT  MODEL 


EHICITA- 


HYPERBOLIC 


OIL  CLEMENT  MODEL 


OIL  ELEMENT  MODEL 
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V.4  Pyke  Cyclic  Simple  Shear 
V.4.1  Motivation 


[Pyke  (1979)]  proposed  a  hyperbolic  model  for  irregular  cyclic  simple 
shear,  which,  unlike  Masing-type  models,  limits  the  peak  shear  stress 
under  arbitrary  loading.  This  is  the  model's  principal  advantage.  Its 
principal  disadvantage  is  its  restriction  to  one-dimensional  simple 
shear.  In  its  present  form,  Pyke's  model  is  not  capable  of  handling 
general  multiaxial  stress-strain  paths,  because: 

a)  The  asymptotic  strength  is  prescribed  as  a  single  magnitude. 

b)  The  tangent  shear  modulus  is  computed  for  simple  shear  about  one 
axis  only. 

c)  There  is  no  provision  for  a  second  tangent  elastic  constant. 

V.4. 2  Assumptions 

The  basic  assumption  of  the  Pyke  cyclic  simple  shear  model  is  that 
the  simple  shear  stress-strain  curve  between  any  two  consecutive  points  of 
shear  strain  reversal  is  a  hyperbola,  with  fixed  initial  slope  immediately 
after  shear  strain  reversal,  and  fixed  upper  and  lower  shear  stress 
asymptotic  limits  of  equal  magnitude. 

V.4. 3  Basic  Equations 

The  general  relation  between  shear  stress  and  shear  strain  for 
irregular  cyclic  simple  shear  is 

_  _  GMAX(y  •  Yc) 


where  y , t 


y»t  =  simple  shear  strain  and  stress 

Yc>  Tc  =  simple  shear  strain  and  stress  at  last  point  of  strain 

reversa 1 

^M/\X  =  slope  of  the  simple  shear  stress-strain  curve 

immediately  after  strain  reversal 

:  y  =  magnitude  of  the  upper  and  lower  shear  stress 

asymptotic  limits 


1  y  =  ;  y  sgn  (d>)  (V.4.2) 

Equations  (V.4.1)  arid  (V.4.2)  are  derived  in  Appendix  0,  and  are  identical 
to  Equations  (0.4)  and  (0.5).  Differentiation  of  Equation  (V.4.1)  yields 
the  tangent  shear  modulus, 


d(T  -  T c ) 


Ty  -  Tc 


[Y  -  Yc) 


(V.4.3) 


Equation  (V.4.3)  is  identical  to  Equation  (0.G). 

V.4.4  Parameter  Determi nation 

The  parameters  and  Ty  are  determined  from  the  virgin  simple 
shear  stress-strain  curve  for  monotonic  loading,  using  a  plot  of  Vy 
versus  i  based  on  the  linear  form  of  Equation  (V.4.1)  when 


which  is 


rc  =  TC  =  °> 


y  =  SlAX  f1  '  T~ 

\  y> 


(V.4.4) 


Equation  (v.4.4)  is  identical  to  Equation  (0.8),  and  the  corresponding 
plot  is  shown  in  Figure  (0.2). 
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Table  (V.4.1)  shows  the  parameters  which  were  used  for  illustrating 
the  behavior  of  Pyke’s  1-D  shear  model.  The  parameters  were  chosen 
somewhat  arbitrarily,  as  no  simple  shear  test  data  is  available  for 
CARES-DRY  Sand. 

V.4.5  Computed  Behavior 

Figure  (V.4.1)  is  the  behavior  predicted  for  simple  shear  by  Pyke's 
model.  The  test  shown  consisted  of  three  fully-reversed  cycles  of  shear 
strain  at  each  of  three  maximum  strain  levels:  0.5  percent,  2  percent, 
and  4  percent.  The  important  features  of  the  model’s  response  are: 

(a)  ratcheting  of  stress- strain  loops  upon  cycling  at  a  given  strain 
1  imi  t 

(b)  softening  of  response  at  progressively  higher  strain  magnitudes 

(c)  stiffened  response  upon  reversal  of  strain  direction 

The  model,  as  it  is  currently  implemented  in  the  SEM,  is  not  general 
enough  to  handle  the  suite  of  exercises  performed  for  the  other  models  in 
this  report. 


TABLE  V.4.1.  PYKE  MODEL  PARAMETERS  FOR  CARES-DRY  SAND 


Parameter 

Symbol 

Variable 

Value 

Units 

Yield  Shear  Stress  y 
Maximum  Shear 

TAUY 

1.0  x  106 

Pa 

Modul us 

fynax 

GMAX 

1.0  x  108 

Pa 

Mass  Density 

p 

PHOREF 

1900 

kg/m^ 

-1 


V.5  El astic-Perfpctly  Plastic 
V.5.]  Motivation 

The  principal  advantages  of  an  elastic-perfectly  plastic  model  are 
that  it  incorporates  stress  state  limits  observed  in  laboratory  strength 
tests,  and  produces  inelastic  strains  when  a  limiting  stress  state  is 
reached.  The  principal  disadvantages  are  that  nonlinear,  inelastic 
behavior  does  not  occur  until  the  failure  surface  is  reached,  and,  with  an 
associated  flow  rule,  predicted  plastic  volume  increases  at  failure  are 
frequently  too  large. 

V.5. 2  Assumptions 

An  elastic-perfectly  plastic,  rate  independent  model  assumes  the 
material  remains  Tinearly  elastic  until  a  limiting  state  of  stress  is 
reached,  defined  by  a  failure  criterion  based  on  strength  test  data. 
Neither  stress  nor  strain  rate  affects  the  stress-strain  relation. 
Inelastic  strains  occur  only  at  failure.  The  flow  rule  can  be  associative 
or  non-associative. 

V.5.?  Basic  Equations 

The  particular  elastic-perfectly  plastic  model  studied  here  has  a 
modified,  associative  Drucker-Prager  failure  criterion  [DiMaggio  and 
Sandler  (1971:94?)],  of  the  form 


ft  -  (*  -  Ce’8Ilj= 


(V.5.1) 


When  the  stress  point  lies  below  the  failure  surface,  so  that 

(V.5. 2) 

the  material  behaves  elastically,  in  accordance  with  the  equations  of 
Appendices  V  and  J.  When  a  stress  increment  calculated  assuming  elastic 
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behavior  moves  the  stress  point  through  the  failure  surface,  so  that 


=  6  >  0 


(V.5.3) 


a  correction  procedure  is  invoked  to  return  the  stress  point  to  the 
failure  surface.  The  correction  procedure  is  developed  in  Appendix  P.  It 


starts  by  computing 


(V. 5.4) 


9Kf  +  Gf 
UTII 


where  K,  G  =  elastic  bulk  and  shear  moduli 


=  -BCe 


-  P  f 1 1 )  E 


(V.5.5) 


(V. 5.6) 


Adjustments  to  the  elastically  computed  invariant  stress  point  coordinates 


are  then  computed  using  the  equations 
^^ADJ  =  (I1}E  "  9KfIdx 


-  Gfjjdx 


I^Vadj 


(V. 5.7) 


(V. 5.8) 


(V.5.9) 


Unless 

f  ( X1  VdJ  » (J^2)ADj]  -  0  (V,J 

the  correction  procedure  is  repeated. 

V.5.4  Parameter  Determi nation 

The  elastic  hulk  and  shear  moduli,  K  and  G,  are  determined  by  the 
methods  described  in  Appendices  V  and  J.  The  failure  criterion 


(V.5.10) 


r.  r.V. . 
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parameters.  A,  B,  and  C,  are  determined  by  nonlinear  regression  or  simply 
by  eye,  using  shear  strength  data  and  the  relations 


1  im 

lj»0 


-BC 


(V.5.11) 

(V.5.1?) 

(V.5.13) 


Parameters  for  the  elastic-perfectly  plastic  model  exercises  are 
shown  in  Table  V.5.1.  The  bulk  and  shear  moduli  and  density  are  equal  to 
those  chosen  for  the  elastic  model.  The  only  difference  between  the 
models  is  the  introduction  of  a  shear  failure  surface. 

V.5.5  Computed  Behavior 

The  behavior  of  the  elastic  and  elastic-plastic  models  is  identical 
for  isotropic  compression  [Figure  (V.5.1)]  and  uniaxial  strain  compression 
[Figures  (V. 5. 1 7- V. 5. 1°) ].  Note  that  the  uniaxial  strain  compression 
stress  path  lies  entirely  under  the  failure  surface  [Figure  (V.5.18)]  for 
these  parameters.  CTC  and  CTE  stress-strain  behavior  [Figures 
( V. 5.2-V. 5.9)3  is  much  improved  over  the  elastic  model  by  the  failure 
surface,  but  only  for  lower  confining  pressures  (o^  <  10  MPa).  A 

non-associated  flow  rule  was  used  here,  so  the  tendency  for  this  material 
to  dilate  with  shearing  is  not  predicted  [Figures  (V. 5.2-V. 5. 5)  for  CTC, 
and  Figures  { V.5.6-V.5.9)  for  CTE], 

The  RTC/E  [Figures  (V.5.10-V.5.13)]  and  PSC/E  [Figures 
(V.5.14-V.5.15)]  calculated  behavior  again  shows  elastic  behavior  until 


the  stress  path  reaches  the  failure  surface.  Behavior  in  the  triaxial 
cell  as  predicted  by  the  elasf ic-perfectly  plastic  model  is  summarized  in 
Figure  (V.5.16). 

Predicted  UXE  behavior  [Figures  (V.5.20)  and  (3.5.21)]  is  more  like 
the  data  than  the  elastic  model  because  of  the  addition  of  the  failure 
surface,  but  still  does  not  match  the  overall  shape  of  the  observed  stress 
path. 

The  results  of  the  axisvmmetric  strain  path  calculations  are  shown  in 
Figures  (V. 5.22)  (WES  paths)  and  ( V.5.23-V.5.24)  (Lade  paths).  In  both 
cases  the  bulk  modulus  is  too  high  in  comparison  with  the  data,  and  so  is 
the  shear  modulus  for  the  Lade  paths.  The  true-triaxial  strain  path  has  a 
similar  problem  as  seen  in  Figures  (V.5.2S)  and  (V.5.26). 


R 


TABLE  V.5.1.  ELASTIC-PERFECTLY  PLASTIC  MODEL 
PARAMETERS  FOR  CARES-DRY  SAND 


Parameter 

Symbol 

Variable 

Val  ue 

Uni  ts 

Bulk  Modulus 

K 

BULK 

3.76  x  108 

Pa 

Shear  Modul us 

G 

SHEAR 

1.44  x  108 

Pa 

A 

CA 

2.88  x  105 

Pa 

Failure  Surface 

B 

CB 

0.00 

1/Pa 

Constants 

C 

CC 

0.00 

Pa 

M 

CAM 

0.215 

Tension  Cu+off 

T 

TCUT1 

6.0  x  10* 

Pa 

Flow  Rule  Switch*  - 

RULE 

1.0 

Mass  Density 

e 

PHOREF 

1900 

kg/m^ 

I 


*0.0  =  Associated 

1.0  =  Non-Associated  (von  Mises  plastic  potential) 
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V.6  Modifier*  AFWL  Engineering 
V.6.1  Motivation 

The  AFWL  engineering  node!  is  hypoel astic-perfectly  plastic  in  shear, 
and  hypoelastic  in  compression.  A  hypoelastic  material  is  one  for  which 
t^e  stress  increments  are  homogenous  linear  functions  of  the  strain 
increments.  The  coefficients  in  the  linear  functions  may  depend  on  stress 
[Fung  (]QPF;44F);  Nelson,  Baron,  and  Sandler  (1071:314)].  The  principal 
advantages  of  the  AFWL  engineering  model  are  ease  of  fitting  to  laboratory 
and  in- situ  test  data,  simplicity  of  the  shear  plasticity  formulation,  and 
the  fact  that  the  model  exhibits  compressive  hysteresis,  which  most  soils 
do  but  many  elastic-perfectly  plastic  models  do  not.  Its  principal 
disadvantages  are  lack  of  hysteresis  in  pure  shear  at  constant  volume 
below  the  failure  surface,  and  lack  of  dilatancy  because  the  plastic 
potential  function  for  yielding  in  shear  is  the  von  Mises  function  (a 
right  circular  cylinder  centered  on  the  hydrostatic  axis),  which  causes 
plastic  incompressibility  in  shear. 

V.6. 2  Assumptions 

Figure  (0.1)  shows  the  shear  failure  surface,  plastic  potential 
surface,  and  hysteretic  hydrostat  which  partially  define  a  modified 
version  of  the  AFWL  engineering  model  in  use  at  Applied  Research 
Associates.  The  initial  loading  hydrostat  governs  only  when  the 
Printout  cancelled  by  operator 


segment  slopes  must  be  specified  so  as  to  avoid  compressive  energy 
creation,  and  Poisson's  ratio  must  be  specified  on  each  hydrostat  segment 
so  as  to  avoid  distortional  energy  creation,  as  well  as  match  constrained 
compression  stress  paths  for  both  initial  loading  and  unloading/ 
reloading.  When  the  compressive  volumetric  strain  falls  below  the  value 
at  which  the  unloading  hydrostat  yields  a  value  of  1^  less  than  T,  the 
tension  cutoff,  both  Ij  and  the  bulk  modulus,  K,  are  set  equal  to  zero. 

As  the  material  recompresses  after  tensile  failure,  both  Ij  and  K  remain 
zero  until  the  compressive  volumetric  strain  exceeds  the  value  at  which 
the  url oadinq/rel oadinq  hydrostat  crosses  the  volumetric  strain  axis.  As 
long  as  K  is  zero  the  shear  modulus,  G,  is  also  zero,  so  that  no 
hypoelastic  deviator  stress  increments  are  generated.  This  means  that  all 
stresses  are  zero,  because  the  only  point  in  Figure  (Q.la)  at  which  the 
tension  limit  can  he  reached  without  previously  violating  the  shear 
failure  surface  is  the  point  where  the  shear  failure  surface  meets  the 
Ij  axis,  and  there  the  deviator  stresses  are  already  zero.  Therefore 
when  tensile  failure  occurs  the  stress  point  automatically  moves  from  the 
tension  cutoff  point  to  the  origin  in  Figure  (Q.la),  and  stays  there  until 
the  hydrostatic  stress  starts  to  build  up  from  zero.  When  that  happens 
the  shear  modulus  again  becomes  nonzero,  and  hypoelastic  deviator  stress 
increments  once  again  start  to  accumulate.  The  AFWL  engineering  model  is 
rate  independent,  and  the  shear  plastic  flow  rule  is  nonassociati ve  with 
respect  to  the  Drucker-  Prage’"  portions  of  the  shear  failure  surface  and 
associative  with  respect  to  the  von  Mises  portion. 
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V.6.B  Basic  Equations 


The  hydrostat  defines  the  incremental  elastic  ( hypoel astic )  bulk 
modulus  as  a  function  of  current  and  maximum  past  compressive  volumetric 
strain,  and  compressive  volumetric  strain  increment. 

K  =  K(eV’cV'1,dEV)  (V.6.1 

where 


=  maximum  past  compressive  volumetric  strain 
Poisson's  ratio  is  also  defined  for  each  hydrostat  segment,  so  that 


v  =  vfrv,£vri,dev)  (V.6.2 

Hypoel astic  constrained  compression  and  shear  moduli  are  then  computer! 
from  the  expressions 


„  3K(1  -  v) 

1  +  v 

,  3K( t  -  ?v) 

G  =  -?cr?  v  ) 


(V.  6.3 

(V  .6.4 


The  shear  failure  surface  is  a  series  of  conical  segments,  each 
having  an  equation  of  the  form 


f  I 


(V  .6.5 


When  the  stress  point  lies  below  the  failure  surface,  so  that 


f  <  0 


(V  .6.6 


the  material  behaves  incrementally  elastically,  in  accordance  with  the 
equations  of  Appendix  J.  When  a  stress  increment  calculated  assuming 
incrementally  elastic  behavior  moves  the  stress  point  through  the  shear 
failure  surface,  so  that 


2  E 


=  5  >  0  (V  .6.7 
a  correction  procedure  is  invoked  to  return  the  stress  point  to  the  shear 
failure  surface.  The  correction  procedure  is  a  special  case  of  that 


developed  in  Appendix  P.  Since  the  von  Mises  plastic  potential  function 


has  the  form 

CQ 

II 

^1 

(  V.6.8) 

it  follows  from  Equations  (V.6.5) 

and  (V.6.R)  that 

(V.6.9) 

fII=  1 

(V.6.10) 

to 

» — i 

II 

o 

(V.6.11) 

9 1 1  "  1 

(V.6.12) 

so  that  Equations  (P.17),  (P.22), 

(P.25),  (P.26) ,  (P.28),  (P.29),  and 

(P.32)  yield 

^3 

dx  = 

(V.6.13) 

M  ■  W 

(V.6.14) 

(dl , )  =  0 

1  P 

(V.6.15) 

(dfi)p  ■  -f3 

(V.6.16) 

(fiy  * a  *  b(ii»3 

(V.6.17) 

(V.6.18) 

Cep  ,  ce  *  M  [sj  {„}T  • 

•§r(s)BT  ,v-6'19) 

From  Equation  (V.6.19),  the  incremental  constrained  modulus  for  a  point  on 
the  shear  failure  surface  is 


n  m  +  3Kb  «. 
I'lw  =  M  +  Si 


(V.6.20) 


and  the  incremental  constrained  horizontal  modulus  for  a  point  on  the 
shear  failure  surface  is 


Mu  =  K  M 
H  o 


3Kb 


37  sls2 


(V.6.21) 


Equations  (V.6.20)  and  (V.6.21)  provide  the  information  required  to 
compute  the  constrained  compression  stress-strain  curve  and  stress  path, 
si  nee 


<ia 

d7 


1 

1 


K 


(V.6.22) 


o 

MV 

=  tt-  (V.6.23) 

0  H 

It  is  shown  in  Appendix  Q  that  the  incremental  stiffness  matrix  defined  hy 
Equation  (V.6.19)  does  produce  a  stress  increment  that  lies  in  the  shear 
failure  surface. 


V.6.4  Parameter  Determination 

The  parameters  of  the  AFWl  engineering  model  are  determined  by 

fitting  a  series  of  straight  lines  to  shear  strength,  hydrostatic 

compression,  and  constrained  (uniaxial)  compression  or  Kq  test  data. 

For  triaxial  compression.  Equations  (K.3)  and  (K.9)  yield 

I,  =  <t  +  2a  (V.6.24) 

l  a  r 

and 


Then  assuming  shear  strength  data  are  obtained  from  drained  triaxial 
compression  tests  at  constant  cell  pressure,  straight  lines  are  fit  to 
consecutive  portions  of  the  data,  plotted  as  (the  dependent 
variable)  versus  ar  (the  independent  variable)  either  by  eye  or  by 
linear  regression.  Successive  linear  relations  between  JjZ  and  I.  are 


then  obtained  by  a  simple  linear  transformation,  and  these  relations 
constitute  the  shear  failure  criterion. 

Since  axial  strain  equals  volumetric  strain  in  a  constrained 
compression  test,  Poisson's  ratio  can  he  calculated  if  both  hydrostatic 
and  constrained  compression  test  data  are  available,  and  if  the  fitted 
straight  line  segments  to  both  curves  have  common  volumetric  strain  break 
points.  Equation  (V.6.3)  yields 


3K  -  M 
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(V.6.26) 


Otherwise,  Poisson's  ratio  is  assumed,  and  either  M  is  calculated  from 
Equation  (V.6.3)  when  only  hydrostatic  test  data  are  available,  or  K  is 


calculated  from  the  equation 

„  M(1  +  v) 

K  = 

when  only  constrained  compression  test  data  are  available. 


(V.6.27) 


In  case  Kq  test  data  are  available  (from  constrained  compression 
tests  conducted  in  a  tri axial  cell,  in  which  the  confining  stress  is 
measured),  Poisson's  ratio  can  be  computed  from  the  hyperbolic  relations 
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(V.6.2P) 


(V.6.29) 


This  last  method  is  the  most  desirable  for  airblast  loading  applications, 
hut  tests  are  more  difficult  and  expensive  than  either  simple 
constrained  compression  (oedometer)  or  hydrostatic  compression  tests. 

Modified  AFV/L  engineering  model  parameters  were  determined  for 
remolded  CAP.ES-DRY  sand  by  fitting  uniaxial  strain  data  (stress-strain  and 
stress  path)  and  triaxial  compression  shear  failure  data.  These 
parameters  are  listed  in  Table  (V.6.1). 
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V.6.5  Computed  Behavior 

By  examining  Figures  (V.6. 17-V.6. 1 9)  one  can  see  that  an  excellent 
fit  is  achieved  for  uniaxial  compression  (UXC),  since  this  data  was  used 
for  fitting.  Isotropic  compression  behavior  (which  is  a  prediction)  is 
also  well  matched,  as  shown  in  Figure  (V.6.1).  CTC  behavior  [Figures 
(V.6.2-V.6.5)]  is  well  fit  with  respect  to  ultimate  shear  strength  [Figure 
(V.6.2)],  but  the  match  of  stress-strain  response  prior  to  failure 
deteriorates  with  increasing  confining  pressures  and  the  model  does  not 
predict  any  dilatancy  [Figures  (V.6.4)  and  (V.6.5)].  Because  the  AFWL 
engineering  model  considers  any  stress  excursion  causing  volume  expansion 
to  be  unloading,  and  subsequently  invokes  a  very  high  stiffness,  CTE 
stress-strain  behavior  is  rather  poorly  predicted  [Figures  (V.6.6)  through 
(V.6.9)].  Note  also  that  a  failure  surface  which  is  symmetrical  about  the 
P-axis  does  not  match  triaxial  extension  data  well.  The  sharp 
discontinuities  in  stress-strain  behavior  caused  by  the  load-unload 
bifurcation  are  again  demonstrated  in  the  RTC/E  calculated  results 
[Figures  (V.6.10-V.6.13)].  Calculated  pure  shear  behavior,  which  by  AFWL 
engineering  standards  is  neither  loading  nor  unloading,  shows  very  stiff 
behavior  [Figures  (V.6.1*)  and  (V. 6. 15)1  because  of  the  Mmgx  convention 
for  deciding  which  bulk  modulus  to  use.  Figure  (V.6.16)  summarizes  the 
behavior  calculated  by  the  AFWL  engineering  model  for  stress  paths  run  in 
the  triaxial  device  starting  at  7.1  MPa  confining  pressure.  An  important 
aspect  of  behavior  to  notice  is  the  single  P-ey  response  enforced  for 
al  1  stress  paths. 

Stress  paths  generated  by  uniaxial  strain  extension  are  not 
characteristic  of  the  data  [Figure  ( V . 6 . 20 ) ]  because  of  the  symmetrical 


failure  surface  and  elastic  response  prior  to  encountering  the  failure 
surface.  Shear  response  for  the  UXE  test  is  shown  in  Figure  (V.6.21). 

Comparisons  of  calculated  strain  path  results  for  this  model  with 
actual  data  are  only  fair.  Stress  path  shapes  for  both  the  WES  strain 
paths  [Figure  (V.6.22)]  and  the  Lade  strain  paths  [Figures  (V.6.23- 
V.6.24)]  are  similar  to  the  data,  but  off  in  both  pressure  and  stress 
difference  magnitudes.  A  large  part  of  the  deviation  occurs  when  the 
specified  strain  path  dictates  expansive  volumetric  strain.  At  this 
point,  as  shown  in  Figure  (V.F.22),  the  shear  modulus  .iumps  due  to  the 
abrupt  bulk  modulus  change  at  unloading.  Eventually  the  failure  surface 
is  reached  and  followed  hack  down  to  its  apex.  The  calculated  response 
for  the  true-triaxial  strain  path  is  compared  with  the  data  from  Nellis 
Baseline  Sand  in  Figures  (V.6.2F)  and  (V.6.2F).  With  the  exception  of  the 
initial  direction  of  the  stress  path,  these  results  are  good. 


TABLE  V.6.1.  MODIFIED  AFWL  ENGINEERING  MODEL 
PAREMETERS  FOR  CARES-DRY  SAND 


Parameter  Symbol  Variable  Value  Uni 


No.  Loaf4  Slope?  n-j 
No.  Unload  Slopes  nu 


Loading  K] 

Bulk  Modul i 


Loading  Strain  e^l 
Break  Points 


Loadi ng 

Poisson's  Ratio  vj 


Unloading 

Bulk  Moduli  K„ 

Unloading 

Hydrostatic  P^u 

Pressure 

Break  Points 

Unloadi ng 

Poi sson' s  Ratio  vu 


Tension  Cutoff  T 

F.S.  Intercept  Y 

F.S.  Slope  S 

Von  Mises  Cutoff  VM 

Mass  Density  p 


RMLS 

8 

- 

RHUS 

5 

_ 

BKL(l) 

0.302  x  107 

P 

BKL(?) 

6.26 1  x  107 

P 

BKL ( 3 ) 

1.491  x  10® 

P 

BKL(4) 

3.530  x  10® 

P 

BKL  ( 5 

1.088  x  109 

P 

BKL(6) 

3.41?  x  109 

P 

BKL  ( 7 

o.04?  x  109 

P 

BKL  f  P ) 

2.102  x  1010 

P 

EBL(l) 

0.01181 

- 

EBL { 2 ) 

0.08191 

- 

EBL ( 3 ) 

0 . 1 2  92 

- 

EBL (4} 

0.164? 

_ 

EBL ( 5 ) 

0.2014 

- 

EPL(6) 

0.2294 

- 

EBL { 7 ) 

0.2516 

- 

EBL ( 8) 

POL ( 1)- 

1.0000 

- 

POL ( 8 ) 

0.32 

. 

BKlKl) 

9.000  x  1011 

P 

BKU(2) 

4.500  x  1010 

P 

BKU(3) 

1.390  x  1010 

P 

BKU  f  4 ) 

4.725  x  109 

P 

BKU(5) 

1.000  x  109 

P 

PBU(l) 

1.045  x  1010 

P 

PBU(2) 

2.800  x  10® 

P 

PBU(3) 

3.000  x  107 

P 

PBU(4) 

2.000  x  107 

P 

PBU(5) 

P0U(1)- 

-2.880  x  10® 

P 

POU { 5 ) 

o.?o 

- 

ST1 

-?.880  x  10® 

P 

Y1 

2.880  x  10® 

P 

SI 

0.215 

- 

VM1 

1.750  x  108 

P 

RHOREF 
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V.7  Effective  Stress  Cap 
V.7.1  Motivation 

The  cap  model,  in  its  various  forms,  sacrifices  some  analytical  and 
computational  simplicity  for  a  more  accurate  representation  of  soil 
behavior  than  provided  by  simpler  models.  The  perfectly  plastic, 
associative  shear  failure  surface  of  the  cap  model  discussed  here 
transitions  exponentially  from  a  Drucker-Prager  asymptote  at  low  confining 
pressure  to  a  von  Mises  asymptote  at  high  confining  pressure.  Volumetric 
hysteresis,  and  control  over  excessive  dilatancy  are  provided  by  a  strain 
hardening,  associative,  ellipsoidal  cap  yield  surface  which  intersects 
both  the  shear  failure  surface  and  the  hydrostatic  axis.  Inside  the 
failure  and  yield  surfaces  the  material  is  hypoelastic  without  hysteresis, 
with  volumetric  and  deviatoric  behavior  uncoupled.  The  cap  model 
satisfies  Drucker's  stability  postulate  and  the  continuity  condition,  and 
thus  produces  stable,  unique  solutions.  The  basic  model  is  a  drained 
(effective  stress)  model,  but  an  undrained  (total  stress)  hydrostat  is 
also  used  to  calculate  pore  pressure  during  undrained  loading.  The 
principal  advantage  of  the  cap  model  is  accuracy  in  representing  most 
aspects  of  soil  stress-strain  behavior.  The  principal  disadvantages  are 
.he  large  number  of  material  parameters  required,  the  amount  of  trial  and 
error  based  on  experience  needed  to  determine  the  parameters,  inability  to 
predict  dilatancy  prior  to  shear  failure,  computational  complexity,  an 
oversimpl ified  approach  to  undrained  response  analysis,  and  lack  of  a 
closed  form  relation  between  total  volumetric  strain  and  effective 
octahedral  normal  stress. 


V.7.2  Assumptions 


The  cap  model  assumes  shear  failure  to  be  governed  by  a  perfectly 
plastic,  associative  failure  surface,  and  plastic  deformation  beneath  the 
shear  failure  surface  to  be  governed  by  a  strain  hardening/softening, 
associative  cap  yield  surface  which  intersects  both  the  shear  failure 
surface  and  the  hydrostatic  axis.  The  cap  strain  hardening  parameter  is 
plastic  volumetric  strain,  which  increases  during  compression  causing  cap 
expansion,  and  decreases  due  to  dilatancy  on  the  shear  failure  surface 
causing  cap  contraction.  Uncoupled  volumetric  and  deviatoric  hypoelastic 
relations  inside  the  failure  and  yield  surfaces  assume  the  bulk  modulus  to 
be  a  function  of  the  octahedral  normal  stress,  and  the  shear  modulus  to  be 
a  function  of  the  octahedral  shear  stress.  Dilatancy  is  recognized  in  the 
drained  (effective  stress)  model  at  failure,  but  not  in  the  undrained 
(total  stress)  model,  which  is  represented  by  an  undrained  hydrostat.  The 
assumption  is  that  total  octahedral  normal  stress  is  uniquely  related  to 
total  volumetric  strain,  regardless  of  shear  deformation.  Pore  pressure 
is  calculated  as  the  difference  between  total  and  effective  octahedral 
normal  stress. 


V. 7.3  Basic  Equations 

The  shear  failure  surface  has  the  equation 

IBI ' ' 

»-e  \ 


=  0 


>  T) 

(V.7.1) 

<  T) 

(V.7.2) 

where  T  is  the  tension  cutoff,  but  the  shear  failure  surface  exists  only 

l  l 

for  Ij  <  k,  where  k  is  the  value  of  1^  at  which  the  shear 

failure  surface  and  the  cap  yield  surface  intersect.  (See  Figure  (R.l).) 


In  the  cap  model  discussed  here  k  is  prevented  from  becoming  negative,  so 


l 

that  k  =  L,  where  L  is  the  value  of  Ij  at  the  center  of  the 
ellipsoidal  cap  yield  surface.  Thus  the  ellipsoidal  cap  has  a  horizontal 
tangent  where  it  meets  the  shear  failure  surface,  and  there  is  no 
von  Mises  transition.  The  cap  yield  surface  equation  is 

J37=  F(I^,k)  =  ^  J(X  -  L)2  -  (I‘  -  L)2  (L<I‘<X)  (V.7.3) 

where 


R 


X  -  L 

=  ~JWT 


(  V. 7.4) 


and  X,  L,  ana  R  are  all  related  to  plastic  volumetric  strain,  as  shown  in 
Figure  (R.4).  In  particular,  X,  the  value  of  1^  at  which  the 
ellipsoidal  cap  yield  surface  intersects  the  hydrostatic  axis,  is  related 
to  plastic  volumetric  strain  by  the  equations 

x  • 3G;  *  i ,n  (jh) 

i 

where  Gr  is  the  hydrostatic  component  of  geostatic  effective  stress,  D 
and  W  are  material  constants,  and 

z  =  fdz  (  V.  7 . 6 ) 

where 

dz  =  0  (k  =  0  and  dej?K  <  o) 

=  (otherwise)  (  V. 7.7) 

Below  the  yield  surfaces  the  material  is  hypoelastic  without  hysteresis, 
with  volumetric  and  deviatoric  behavior  uncoupled.  Thus 
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( V.7.10) 


(V.7.11) 


and  Kis>  Kis*  k2s’  Gi*  Gl*  and  G2  are  material  constants. 

Equations  (V.7.1)  through  (V.7.11)  apply  to  the  drained  (effective 

stress)  cap  mode).  In  addition,  an  undrained  (total  stress)  bulk  modulus 
for  isotropic  compressive  loading  is  defined  by  an  equation  identical  in 
form  to  Equations  (V.7.10)  and  (V.7.11). 
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In  an  undrained  calculation  pore  pressure  is  calculated  by  the  equation 
1 1  -  i; 

u  =  >-3  1  ( V.  7.13) 


A  more  detailed  mathematical  description  of  the  cap  model  is  given  in 
Appendices  R  and  S. 

V.7.4  Parametric  Determination 

All  the  cap  model  parameters  can  in  principle  be  determined  from 
standard  laboratory  isotropic  compression  and  triaxial  compression  tests 
at  constant  cell  pressure,  as  explained  in  Appendicies  R  and  S.  However, 
the  parameters  thus  determined  are  usually  refined  by  trial  and  error  to 
match  uniaxial  and  triaxial  compression  test  data,  and  sometimes  even 
dynamic  field  test  data.  The  trial  and  error  process  is  based  on 
experience,  and  has  not  been  explained  in  step  by  step  fashion  in  the 
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1 iterature. 


Table  (V.7.1)  lists  the  cap  model  parameters  used  for  modeling 
remolded  CARES-DRY  Sand.  The  functions  for  the  elastic  bulk  and  shear 
moduli  have  been  reduced  to  constant  linear  relationships  (K..,  ) . 

This  is  because  when  K  =  K ( I ^ )  and  G  =  G  1  ^  1S  difficult  to 

maintain  a  reasonable  uniaxial  strain  unload-reload  stress  path  shape. 

The  failure  surface  parameters  (Cg  a)  were  fit  to  standard  triaxial  test 
results  on  dry  material.  The  cap  shape  (R)  was  chosen  based  on  previous 
experience  with  dry  alluvium,  and  its  hardening  parameters  (W,  D)  were 
iteratively  fit  to  uniaxial  compression  stress-strai n  data. 

V.7.5  Computed  Behavior 

V.7.5.1  Drained  (or  Dry)  Model 

As  noted  above,  uniaxial  compression  stress-strain  data  was  used  for 
fitting  the  cap  model.  Figure  (V.7.18)  shows  this  data,  compared  with 
computed  behavior.  The  model  cannot  match  the  changes  in  curvature  of  the 
stress-strain  data  because  of  the  single  exponential  formulation  of  the 
cap  hardening  function  [Equation  (R.?l)].  Rather,  it  is  fit  to  produce  an 
acceptable  response  over  the  stress  range  of  interest  (ofl  -  o  -  50  MPa, 
in  this  case).  Stress  path  data  from  the  Kq  test  [Figure  (V.7.19)]  is 
matched  reasonably  well,  but  this  model  does  not  produce  any  unload-reload 
loops  due  to  the  purely  elastic  behavior  under  both  the  failure  surface 
and  cap.  The  variation  of  shear  stiffness  in  the  uniaxial  test  is  shown 
in  Figure  (V.7.?0).  A  comparison  of  the  calculated  isotropic  compression 
response  with  test  data  [Figure  (V.7.1)]  is  similar  to  the  calculated  vs. 
observed  DXC  comparison.  Overall,  the  fit  is  good,  but  sharp  changes  in 
stiffness  with  increasing  pressure  cannot  be  matched. 


Calculated  CTC  stress-strain  response  is  better  at  lower  confining 
pressures  than  at  the  higher  ones  [see  Figures  (V.7.2)  and  (V.7.3)].  The 
tests  at  59  MPa  and  100  MPa  confining  pressure  show  much  too  stiff  a 
response  but  with  a  reasonable  stress  difference  at  failure.  This  is 
because  at  higher  initial  confining  pressures,  the  capacity  for 
irrecoverable  volumetric  compaction  has  nearly  been  exhausted,  causing 
very  stiff  volumetric  response  [see  Figure  (V.7.5)].  Figure  (V.7.4)  shows 
that  no  dilation  is  predicted  by  the  cap  model  for  CTC  because  the  stress 
point  is  alv'ays  located  on  the  cap  or  at  the  intersection  of  the  cap  and 
failure  surface.  Here,  the  plastic  strain  rate  vector  is  forced  to  be 
perpendicular  to  the  P-axis. 

CTF  response  [Figures  (V. 7.6-V. 7. 9) ]  at  low  strains  is  very  stiff, 
because  initially  the  cap  does  not  need  to  expand  in  order  for  the  stress 
point  to  move  away  from  the  P-axis.  Thus,  this  initial  behavior  is 
elastic  [Figures  (V.7.6)  and  (V.7.7)].  However,  further  outward  movement 
of  the  stress  point  requires  the  cap  to  move  out  to  maintain  its  shape 
(R  =  2.5).  This  cap  hardening  is  accompanied  by  compressive  volumetric 
strain  [Figure  (V.7.9)]  and  a  much  softer  stress-strain  response. 
Eventually,  the  shear  failure  surface  is  reached,  whereupon  volumetric 
straining  stops  [see  Figure  (V.7.8)].  Thus  there  is  no  predicted  dilation 
here  either,  and  in  fact,  the  data  tends  to  support  this  [Figure  (V.7.8)]. 

The  RTC  and  RTE  exercises  yield  very  different  behavior  from  each 
other,  as  shown  in  Figures  (V.7,10)  through  (V.7.13).  During  the  RTC 
test,  the  stress  point  moves  off  the  P-axis  into  the  elastic  region  and 
thus  yields  very  stiff  Initial  stress-strain  behavior.  The  failure 
surface  is  soon  encountered,  however,  and  the  material  fails  in  shear  and 


begins  to  dilate  [Figure  (V.7.12)]  at  constant  pressure  [Figure  (V.7.13)]. 
The  cap  subsequently  retracts  to  meet  the  current  stress  point  and  dilation 
is  then  stopped.  RTE  behavior,  conversely,  is  initially  soft  because  the 
cap  is  being  forced  outward  constantly  [Figure  (V.7.10)].  RTE  volume 
response  is  purely  compressive.  In  fact,  at  15  percent  radial  strain  in 
this  test,  the  shear  failure  surface  has  still  not  been  reached  at  any  of 
the  confining  pressures. 

The  PSC  and  PSE  stress  paths  [Figures  (V.7.14-V.7.16)]  both  force  the 
stress  point  directly  outward  from  the  P-axis,  causing  cap  expansion  and 
therefore  volumetric  compression  [Figure  (V.7.16)].  When  the  failure 
surface  is  reached,  the  plastic  strain  rate  vector  swings  perpendicular  to 
the  P-axis  and  volumetric  strain  ceases. 

Figure  (V.7.17)  summarizes  the  behavior  of  the  cap  model  for  the  tests 
performed  in  the  triaxial  cell  starting  at  7.1  MPa  confining  pressure.  Cap 
model  uniaxial  extension  (UXE)  behavior  is  elastic  until  shear  failure, 
producing  a  stress  path  substantially  different  from  that  actually  observed 
[Figure  (V.7.21)].  Figure  ( V. 7.22)  shows  that  dilation  occurs  when  the 
stress  point  is  on  the  failure  surface,  meaning  that  the  cap  retraction  in 
this  case  is  not  fast  enough  to  overtake  the  stress  point  and  limit 
subsequent  dilation. 

Cap  model  predictions  of  axisymmetric  strain  path  experiments  are 
shown  in  Figures  ( V. 7.23)  through  (V.7.25).  The  initial  behavior  in  these 
tests  is  dominated  by  an  adjustment  from  isotropic  stress  to  uniaxial 
strain  conditions.  Apparently,  the  stress  point  moves  up  the  cap  until  a 
point  on  the  cap  is  reached  which  has  an  outward-normal  plastic  strain-rate 
vector  which  produces  a  compressive  increase  in  radial  stress  under 


uniaxial  strain  boundary  conditions.  Until  this  point  is  reached,  radial 
stress  increments  are  tensile  and  the  mean  normal  stress  drops.  When  this 
point  is  reached,  the  stress  path  assumes  a  slope  more  typical  of  uniaxial 
total  strain.  From  there,  the  stress  path  breaks  over  toward  the  failure 
surface  when  the  strain  path  dictates  volume  expansion,  and  subsecjuently 
follows  the  failure  surface.  With  the  exception  of  this  initial  behavior 
the  cap  model  matches  both  the  WES  data  [Figure  (V.7.23)]  and  the  Lade 
data  [Figures  (V.7.24)  and  (V.7.25)]  reasonably  well. 

Figures  (V.7.26)  and  (V.7.27)  compare  the  truly-triaxial  strain  path 
prediction  for  remolded  CARES-DRY  Sand  with  data  from  Ko's  tests  on  Nellis 
Baseline  Sand  [Ko  and  Meier  (1983)].  The  major  departures  from  typical 
alluvium  data  occur  again  at  initial  departure  from  the  P-axis  and  during 
reloading  from  the  spalled  condition. 

V.7.5.2  Undrained  Model 

Since  no  data  exists  for  undrained  tests  on  saturated  remolded 
CARES-DRY  Sand,  only  calculated  results  for  CTC/CTE  will  be  shown.  This 
will  serve  to  illustrate  how  the  undrained  portion  of  the  effective  stress 
cap  model  works. 

Figure  (V.7.28)  shows  the  total  and  effective  stress  paths  for  the 
CTC  and  CTE  tests  at  three  confining  pressures  (7,  59,  and  100  MPa).  The 
difference  between  the  two  paths  is  the  predicted  pore  pressure,  which  is 
plotted  against  axial  strain  in  Figure  (V.7.29).  Stress-strain  response 
for  the  axial  and  radial  directions  is  shown  in  Figure  (V.7.30)  and  the 
volume  strain  prediction  is  plotted  against  axial  strain  in  Figure 
(V.7.31).  Note  that  the  volume  strain  in  Figure  (V.7.31),  while  not  zero 
for  undrained  loading,  is  almost  two  orders  of  magnitude  less  than  that  in 
Figure  (V.7.1). 


TABLE  V.7.1.  EFFECTIVE  STRESS  CAP  MODEL  PARAMETERS 
FOR  REMOLDED  CARES-DRY  SAND 


Parameter 

Symbol 

Variable 

Val  ue 

Units 

Drained  Elastic 

Ki 

AKI 

*.0  x  109 

Pa 

Bulk  Modulus 

*1 

AK1 

-0- 

*2 

AK2 

-0- 

Undrai  ned 

Kim 

AKIM 

8.0  x  109 

Pa 

Bulk  Modulus 

Kin 

AKIM 

-0- 

K2n 

AK2M 

-0- 

Gi 

AGI 

3.0  x  109 

Pa 

Shear  Modul us 

G1 

AG1 

-0- 

g2 

AG? 

-0- 

Cs 

AC 

2.88  x  105 

Pa 

Failure  Surface 

a 

AM 

0.215 

-- 

B 

RB 

-0- 

C 

CCC 

-0- 

Ri 

ARI 

2.50 

^  _ 

Cap  Shape 

Ri 

AR1 

-0- 

R2 

AR2 

-0- 

Cap  Hardening 

w 

AW 

0.200 

_  _ 

D 

AD 

1.800  x  io-fi 

1/Pa 

Mass  Density 

p 

RHOREF 

1900 

kg/n^ 
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V.  8  Lade 


V.8.1  Motivation 

The  Lade  model  discussed  here  is  an  el astopl astic  model  with  two 
yield  surfaces.  One,  called  the  expansive  yield  surface,  is  bullet  shaped 
with  its  nose  at  the  origin  in  stress  space.  The  other,  called  the 
collapse  yield  surface,  is  spherical  with  its  center  at  the  origin.  Both 
yield  surfaces  harden  in  response  to  the  corresponding  plastic  work,  and 
the  expansive  yield  surface  also  softens  when  the  corresponding  plastic 
work  exceeds  a  certain  value.  The  collapse  yield  surface  is  associative 
and  the  expansive  yield  surface  non-associative.  The  principal  advantage 
of  the  Lade  model  is  accuracy  in  representing  most  aspects  of  soil  stress- 
strain  behavior.  The  model  exhibits  nonlinear,  inelastic  behavior  in  both 
shear  and  compression  even  at  small  strains,  and  the  expansive  yield 
surface  has  a  non-circular  octahedral  cross-section  and  therefore 
indicates  an  influence  of  the  intermediate  principal  effective  stress  (or 
Lode's  parameter)  on  shear  strength.  The  principal  disadvantages  of  the 
Lade  model  are  possible  underprediction  of  compressibility  under  the 
influence  of  shear  at  small  strains,  lack  of  flexibility  in  matching  true 
triaxia1  shear  strength  data  in  the  octahedral  plane,  lack  of  a  device  to 
prevent  negative  plastic  work,  and  possible  instability  and  lack  of 
uniaueness  due  to  strain  softening  of  the  expansive  yield  surface. 

V.8.2  Assumptions 

The  assumptions  underlying  the  Lade  model  are  basically  those 
discussed  in  Appendix  0,  except  the  dissipation  condition  is  not 
enforced.  Instead  of  a  yield  function  defining  the  stress  levels  at  which 
plastic  strain  increments  can  occur,  it  defines  the  stress  levels  at  which 
plastic  strain  increments  will  occur  [Lade  and  Nelson  (1981 rSO*, 507)]. 
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The  strain  hardening  parameter  for  each  yield  surface  is  the  corresponding 
plastic  work,  rather  than  plastic  volumetric  strain  as  with  the  cap 
model-  The  use  of  two  separate  yield  surfaces,  each  with  its  own 
potential  surface  and  hardening  rule,  assumes  that  a  material  really  has 
two  separate  yield  surfaces. 

V.8.3  Basic  Equations 

The  collapse  yield  criterion  has  the  equation 

fc  =  -  f‘‘ .  o  (v.8.1) 


i n  whi ch 


fc  '  *1  212  '  3(o0CT  +  1  OCT* 


f"  n2(WcV/P 


(V.8.2! 


(V.8.3] 


where 


Pa  =  atmospheric  pressure 
C,  p  =  material  parameters 


“c  *  /htm 


(V  .8.4) 


The  collapse  plastic  potential,  to  which  the  collapse  plastic  strain 


increment  vector,  jdec]  ,  is  normal,  is 
Sc  *  fc  -  >1  *  2I2 


(  V.8.51 


The  expansive  yield  criterion  has  the  equation 


Wi 


i n  whi ch 


where 

m»  ^2 »  a,  b,  q  =  material  parameters 
and 

-p  ■  flV  M 


(  V.8.7) 


(  V.8.8) 


(  V. 8.9) 


(  V.8.10) 


The  expansive  plastic  potential,  to  which  the  expansive  plastic  strain 
increment  vector,  [d£Pl ,  is  normal,  is 


9 


P 


(  V.8.11) 


where  n2  is  also  a  material  parameter. 

The  uni oadi ng/rel oadi ng  elastic  modulus  is  assumed  to  be  given  by  the 
expression 


E  =  K  p  f—)n 
ur  urHa  \  Pd  / 


(V  .8.12) 


where 


Kur,  n  r  material  parameters 
Poisson's  ratio  is  usually  assumed  constant. 

A  method  for  calculating  the  octahedral  cross-section  of  Lade's 
failure  surface  (the  expansive  yield  surface  at  its  maximum  extent)  is 
given  in  Appendix  T. 


The  computational  features  of  Lade's  model  are  covered  in  Appendices 
F,  G,  H,  I,  J,  K,  and  T.  Using  the  value  of  Young's  modulus  from  Equation 
(V.8.12)  and  an  assumed  Poisson's  ratio,  the  elastic  incremental  stiffness 
matrix  is  given  by  Equation  (J.31).  The  general  equations  of 
elastoplasticity  given  in  Appendix  D  apply,  except  that  the  dissipation 
condition,  Equation  (D.8)  is  not  enforced.  The  el astopl asti c  incremental 
stiffness  matrices  are  calculated  as  described  in  Appendix  G,  but  the 
polar  mode  check  described  in  Appendix  I  is  not  used.  Instead,  if  the 
stress  point  lies  on  a  yield  surface,  that  surface  is  assumed  to  be 
active.  However,  only  positive  plastic  work  increments  are  accumulated  in 
calculating  total  plastic  work. 

V.8.4  Parameter  Determination 

Determination  of  the  material  parameters  in  Equations  (V.8.1)  through 
(V.8.12)  by  a  series  of  linear  laboratory  test  data  plots  is  explained  in 
Appendix  T.  Parameters  for  CARES-DRY  sand  were  determined  by  using  an 
automatic  fitting  routine  in  the  Soil  Element  Model  which  is  based  on 
Appendix  T.  Step-by-step  results  from  this  fitting  process  are  shown  in 
Figures  (V.8.1)  through  (V.8.18). 

Elastic  stiffness  data  (Young's  Modulus)  as  determined  from  the 
unload-reload  portions  of  isotropic  compression,  uniaxial  strain,  and 
triaxial  compression  tests  presented  by  [Cargile  (1984)]  are  plotted 
versus  radial  confining  pressure  in  Figure  (V.8.1).  This  plot  yields 
values  for  the  modulus  number,  «ur,  and  exponent,  n. 

Figure  (V. 8.2)  shows  the  loading  hydrostatic  data  used  for  obtaining 
collapse  yield  surface  parameters.  In  Figure  (V.8.3)  the  hydrostat  has 
been  integrated  to  show  collapse  plastic  work  versus  radial  confining 
pressure  squared.  Elastic  properties  are  varied  according  to  the 


previously  determi ned  relationship  with  during  integration.  A  single 
straight  line  is  passed  through  the  data  (open  circles)  to  determine  the 
collapse  constant,  C,  and  the  exponent,  p. 

Figure  (V.8.4)  shows  triaxial  stress-strain  data  at  five  confining 
pressures.  When  the  data  does  not  approach  a  "peak"  response,  it  needs  to 
be  artificially  extended  by  some  method.  Extension  to  peak  is  necessary 
for  this  type  of  hardening  model  to  produce  reasonably  shaped 
stress-strain  curves  having  the  correct  value  of  expansive  plastic  work  at 
peak  stress  difference.  The  method  arbitrarily  chosen  here  was  to  assume 
a  hyperbolic  shape  for  the  last  twenty-five  percent  of  both  the  axial  and 
radial  stress  di fference-strain  curves.  A  Kondner  plot  of  strain  divided 
by  stress  difference  vs.  strain  is  then  used  to  extend  the  data  to 
ninety-five  percent  of  peak  stress  difference  as  predicted  by  the  straight 
line  fit.  Two  such  plots  are  shown  in  Figures  (V.8.5)  and  (V.8.6)  for  the 
CARES-DRY  data  at  3.4  MPa  and  7.0  MPa,  respectively. 

Peak  response  data  from  the  extended  stress-strain  curves  are  plotted 
in  Figure  (V.8.7)  as  fp>max  vs.  pa/ I ^ .  The  parameters  from  this 
fit,  nj  and  m,  define  the  shape  and  most  expanded  position  of  the 

l 

expansive  yield  surface,  f 

Given  the  previously  determined  elastic  and  collapse  plastic  work 
relationships,  collapse  and  expansive  plastic  strains  resulting  from  the 
triaxial  test  are  calculated  and  plotted  in  Figures  (V.8.8)  and  (V.8.9), 
respecti vely.  Figures  (V.8.10)  through  (V.8.14)  show  the  components  of 
volume  strain  which  are  relevant  to  the  Lade  model:  total,  elastic, 
collapse  plastic,  and  expansive  plastic  (which  is  the  difference  between 
total  and  elastic  +  collapse).  The  basic  theory  of  the  model  requires 
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that  the  strain  associated  with  the  expansive  yield  surface  always  be 
expansive.  Rut  in  Figures  (V.R.ll),  (V.P.l?),  and  (V.8.13),  this  does  not 
always  happen.  The  reason  for  this  is  that  at  these  confining  pressures, 
collapse  behavior  as  determined  by  the  isotropic  test  is  stiffer  than  would 
be  indicated  by  the  triaxial  test.  One  remedy  for  this  inconsistency  is  to 
use  an  envelope  of  pressure-vol ume  response  which  encompasses  both  tests. 

The  triaxial  data  and  calculated  expansive  plastic  strains  are  used  to 

l 

determine  the  variation  of  the  expansive  yield  surface,  f  ,  with 
expansive  plastic  work,  W  This  is  shown  in  Figure  (V.8.15).  Note  that 
the  data  does  not  conform  to  the  idealized  shapes  for  these  curves  as 
postulated  by  the  model.  Values  of  plastic  work  at  peak  stress  difference 
lay  roughly  on  a  straight  line  in  log  space  [Figure  (V.8.16)],  but  the  shape 
parameter,  a,  in  arithmetic  space  does  not  [Figure  (V.8.17)]. 

The  expansive  plastic  potential  parameters  are  derived  from  Figure 

II 

(V.8.1P),  which  shows  the  variation  of  with  both  f  and  o^. 

The  straight  lines  lie  in  a  plane  which  passes  through  the  entire  set  of 
data.  The  rough  appearance  of  the  data  is  primarily  a  result  of  three 
factors: 

(a)  The  data  was  digitized  by  hand  from  small  plots,  and  radial  strain 
was  calculated  by  graphically  subtracting  strain  difference  from 
axial  strain. 

(b)  (Inload-reload  loops  were  eliminated  from  the  data,  but  did  in  fact 
have  a  significant  effect  on  volumetric  response. 

(c)  the  fitting  process  up  to  this  point  may  introduce  some 
inconsistency  between  test  data  and  predicted  model  response. 

If  raw  data  had  been  available  in  digital  form,  and  only  loading  had  been 
performed  in  these  specific  tests,  a  much  smoother  plot  would  be  expected. 


Table  (V.8.1)  summarizes  the  Lade  model  parameters  for  remolded 


CARES-DRY  sand. 

V.8.5  Computed  Behavior 

With  a  two-parameter  exponential  fit  to  the  hydrostat,  the  Lade  model 
can  produce  a  curve  of  single  inflection  only,  either  convex  or  concave. 
Figure  (V.8.19)  shows  this,  compared  with  the  test  data.  The  fit  is 
reasonable  at  higher  stresses,  but  cannot  match  the  low  stress  variations 
in  bulk  modulus  trend. 

CTC  stress- strain  data,  used  to  fit  the  model,  is  well  matched 
[Figures  (V.8.20)  and  (V.8.21)].  Volumetric  response  is  qualitatively 
good  [Figure  IV. 8. 22)],  but  overestimates  the  tendency  for  dilation  in  the 
looser  (lower  a^)  sand?.  Pressure-vol ume  response  comparisons  [Figure 
(V.8.23)]  are  very  good,  again  with  the  exception  of  over-dilation  at  low 
stresses.  CTE  shear  failure  levels  [Figure  (V.8.24)]  could  be  matched 
better  with  a  slightly  lower  asymmetry  in  the  expansive  yield  surface. 
Initial  shear  stiffness  is  too  high  [Figure  (V.8.2F)].  The  test  specimens 
tended  to  compact  only,  while  the  model  predicts  a  small  amount  of  initial 
elastic  expansion,  subsequent  plastic  compaction  as  the  collapse  yield 
surface  is  pushed  out,  and  finally  dilation  as  the  expansive  yield  surface 
is  pushed  toward  its  ultimate  position  [Figures  (V.8.26)  and  (V.8.27)]. 

The  constant  axial  stress  (RTC/RTE)  and  constant  mean  normal  stress 
(PSC/PSE)  tri axial  test  exercises  serve  to  illustrate  the  principal 
aspects  of  the  Lade  model:  work  hardening/softening,  asymmetric  yield 
surface,  and  shear-volume  coupling  (dilatancy).  Due  to  the  work-hardening 
expansive  yield  surface,  stress-strain  behavior  is  smooth  for  both  the 
RTC/E  tests  [Figures  (V.8. 2P-V.8. 31 ) ]  and  PSC/E  tests  [Figures  ( V . 8 . 32- 
V.8.34)].  Different  levels  of  stress  difference  near  shear  failure  for 
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PSC  and  PSE  [Figure  (V.8.32)]  are  a  result  of  the  asymmetry  of  the  failure 
surface.  Work  softening  is  most  evident  in  the  RTE  test  [Figure  (V.8.28)]. 
Note  that  this  causes  a  negative  net  shear  stiffness,  as  seen  in  Figure 
( V.8.29) .  Dilation  is  evident  for  both  the  RTC/RTE  tests  [Figures  (V.8.30) 
and  (V. 8. 31)1  and  the  PSC/PSE  tests  [Figure  (V.8.34)].  This  is  especially 
important  in  the  constant  pressure  tests. 

Behavior  predicted  by  the  Lade  model  for  tests  run  in  the  triaxial 
device  from  7.1  MPa  confining  pressure  is  summarized  in  Figure  (V.8.35). 
Note  the  variability  in  pressure  volume  response  due  to  varying  dilation. 

Predicted  uniaxial  strain  compression  (UXC)  behavior  compares  poorly 
with  thp  data  [Figures  (V.8.38)  through  (V.8.38)].  The  calculated  uniax  is 
very  soft  at  low  axial  stress,  due  to  the  low  initial  confining  pressure 
(0.1  MPa)  and  subsequent  low  elastic  stiffness  (Eur  is  a  function  of 
radial  stress).  The  model  also  suddenly  softens  at  about  os  =  38  MPa 

d 

[Figures  (V.8.36)  and  (V.8.37)]  because  the  peak  expansive  plastic  work  is 
attained.  Clearly,  this  type  of  softening  behavior  is  not  appropriate  for 
this  material  under  these  test  conditions.  Qua! itatively,  the  UXE  stress 
paths  compare  onite  well  with  thp  data  [Figure  (V.8.39)],  but  again  it 
appears  that  somewhat  more  shear  capacity  in  extension  is  required.  Figure 
(V.8.40)  shows  calculated  UXE  stress- strain  response. 

Calculated  results  for  the  WES  strain  paths,  shown  in  Figure  (V.8.41) 
are  reasonable  hut  show  a  slightly  high  shear  stiffness,  as  can  be  seen  in 
the  plot  of  ( n ^ - r ^ )  vs.  (ej-ej).  Neither  calculated  stress  path 
shows  the  rather  flat  stress  difference  response  upon  increased  radial 
expansion  as  is  observed  in  the  data.  The  second  set  of  strain  paths 
(denoted  "Lade")  originates  at  a  low  confining  pressure  (o3c  =  0.4  MPa). 
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This  causes  the  Lade  model  to  have  a  very  low  initial  volumetric  stiffness 
and  causes  underestimated  pressures  [Figures  (V.8.42)  and  (V.8.43)].  The 
overall  shape  of  the  calculated  stress  paths  is  very  good.  Figures 
(V.B.44)  and  (V.8.4F)  show  predicted  response  for  the  true  triaxial  strain 
path.  Data  for  Nellis  Baseline  Sand  is  included  on  the  plots  for 
comparison.  Both  the  stress  path,  and  volume  response  comparisons  look 
quite  reasonable. 
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TABLE  V. 8. 1(a).  LADE  MODEL  PARAMETERS  FOP 
REMOLDED  CARES-DRY  SAND 


Parameter  Symbol  Variable  Value  Units 


Modulus  Number 

K 

u  r 

EKl'P 

383.5 

-- 

Modulus  Exponent 

n 

EN 

0.8412 

Poisson's  Patio 

V 

POIS 

0.20 

-- 

Collapse  Constant 

r 

C 

2.70 0  x  10‘3 

-- 

Collapse  Exponent 

p 

PC 

0.6482 

-- 

Yi el d  Constant 

nl 

ETA1 

84.5? 

-- 

Yield  Exponent 

m 

CURVM 

0.2261 

-- 

r 

R 

0.1806 

Plastic  Potential 

s 

ss 

0.7309 

-- 

Constants 

t 

T 

-14.52 

-- 

a 

ALPHA 

1.Q72 

Work-Hardening 

e 

BETA 

2.140  x  10‘4 

-- 

Constants 

P 

PV! 

0.6379 

-- 

1 

ELW 

0.8077 

Peak  Plastic 

P 

Wp,pk 

0 

WPPK 

Pa 

Exp.  Work 

a 

A 

varies,  see 

__ 

b 

B 

Table  V . 8. 1  (b) 

-- 

Initial  Confining 

n3c 

SIGMA 3 

MPa 

Pressure 

Initial  Plastic 

w  n 

WC 

Pa 

Comp.  Work 

Mass  Density 

C  j  0 

p 

RHOREF 

1900 

kg/m' 
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TABLE  V. 8. 1(b).  LADE  MODEL  PARAMETERS  WHICH  VARY  WITH  CONFINING  PRESSURE 


°3c 

Q 

wppk 

A 

B 

WC,0 

(MPa) 

(Pa) 

(Pa) 

0.1 

1.97? 

6.40 

X 

104 

177.3 

7.928 

X 

10'6 

5.48  x 

102 

0.4 

1.97? 

2.0? 

X 

105 

99.06 

2.516 

X 

10"6 

3.46  x 

103 

1.8 

1.976 

6.60 

X 

105 

54.31 

7.664 

X 

10~7 

2.3?  x 

104 

3.5 

1.97° 

1.13 

X 

106 

41.44 

4.471 

X 

10'7 

5.50  x 

104 

7.0 

1.987 

1.98 

X 

lO6 

31.34 

2.545 

X 

10~7 

1.35  x 

105 

20.0 

2.014 

4.6? 

X 

106 

20.85 

1.075 

X 

10~7 

5.27  x 

105 

32.0 

2.040 

6.75 

X 

106 

17.6? 

7.264 

X 

10'8 

9.70  x 

105 

59.0 

2. 0°6 

1.10 

X 

107 

14.54 

4.324 

X 

lO'8 

2.13  x 

106 

100.0 

2.183 

1.69 

X 

107 

12.81 

2.704 

X 

ID'8 

4.25  x 

106 

S  I 

[  I 


I 
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FIGURE  V.8,20  LADE  MODEL.  EXERCISE-TRIAXIAL  COMP(CTC)  -  STRESS  DIFF  VS.  AXIAL  STRAIN 


AEOER  SOIL  ELEMENT  MODEL 


VOLUMETRIC  STRAIN 


OIL.  ELEMENT  MODEL 


;OIL  CL  EMENT  MODEL 


Q-i  OS  S 'Z  O’O  VZ~  O’S- 

pi*  (vd>  jjia  ssBdis 


VL- 


O'Ol-  S'Cl- 


531 


EAXIAL  &  ERADIAL 

FIGURE  V. 8.28  LADE  MODEL  EXER-REDUCED  TRIAX(RTC+RTE)  -  STRESS  DIFF  VS.  STRAIN 
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AXIAL  STRAIN 

FIGURE  V.8.30  LADE  MODEL  EXER-REDUCED  TRIAX(RTC+RTE)  -  AXIAL  STRAIN  VS  VOLUME  STRAIN 
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FIGURE  V, 8, 31  LADE  MODEL  EXER-REDUCED  TRIAX(RTC+RTE)  -  PRESSURE  VS.  VOLUMETRIC  STRAIN 
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AXIAL  STRAIN 

FIGURE  V . 8. 34  LADE  MODEL  EXERCISE-PURE  SHEAR(PSC+PSE)  -  AXIAL  STRAIN  VS  VOLUME  STRAIN 
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APPENDIX  W 


DEVELOPMENT  OF  THE  ARA  CONIC  MODEL 

W.l  Motivation 

Of  the  eight  soil  constitutive  models  examined  in  Section  3  and 
Appendix  V,  the  Lade  model  is  most  appealing  from  three  important 
standpoi nts: 

a)  favorable  rating  with  respect  to  seven  of  the  ten  evaluation 
criteria  in  Table  3.1; 

b)  accuracy  and  flexibility  in  representing  soil  stress-strain 
behavior;  and 

c)  ease  of  developing  intuition  for  parameter  physical  significance 
and  accuracy. 

Consequently  ARA  elected  to  modify  the  Lade  model  rather  than  create  a 
completely  new  one,  to  develop  a  soil  constitutive  model  suitable  for 
analyzing  the  response  of  soil  masses  to  complex  dynamic  loadings. 

The  modifications  were  designed  to  achieve  the  following  additional 
desirable  features: 

a)  better  volumetric  strain  response  under  non-i sotropic  loading; 

b)  greater  flexibility  in  matching  shear  strength  data,  in  both  the 
triaxial  and  octahedral  planes; 

c)  correct  plastic  mode  selection  based  on  the  thermodynamically 
related  dissipation  condition  that  a  positive  plastic  work 
increment  accompany  yielding; 

d)  finite,  reasonable  friction  angle  at  low  confining  pressure; 

e)  essentially  constant  shear  strength  at  high  confining  pressure; 
and 

f)  direct  (noniterative)  shear  strength  calculation  in  both  the 


triaxial  and  octahedral  planes. 

Several  Lade  model  features  have  been  retained: 

a)  the  basic  model  construction,  i.e.,  two  yield  surfaces,  one 
compressive  and  one  expansive,  both  strain  hardening,  the 
compressive  yield  surface  associative  and  the  expansive  yield 
surface  non-associetive; 

b)  both  the  compressive  and  expansive  work  hardening  formulations; 
and 

c)  the  unloading/reloading  elastic  modulus  formulation. 

New  features  include: 

a)  an  ellipsoidal  compressive  yield  surface  to  increase 
compressibility  in  the  presence  of  shear  deformation; 

b)  a  hyperbolic  expansive  yield  surface  vnth  a  triple  ellipsoidal 
octahedral  cross  section,  possessing  a  finite,  adjustable  slope 
(friction  angle)  at  low  confining  pressure,  essentially  constant 
shear  strength  at  high  confining  pressure,  flexibility  in 
matching  both  compression  and  extension  shear  strength  data,  a 
completely  smooth  octahedral  cross  section,  and  directly 
computable  shear  strength; 

c)  enforcement  of  the  dissipation  condition; 

d)  development  of  a  polar  mode  check  based  on  the  dissipation 
condition,  to  determine  uniquely  and  without  trial  and  error 
which  yield  surfaces  are  active  under  a  given  state  of  stress  and 
prescribed  total  strain  increment;  and 

e)  determination  of  compressive  yield  surface  parameters  by  fitting 
the  plastic  hydrostat  directly  (using  a  linear  transformation) 
rather  than  having  to  compute  compressive  plastic  work. 


In  addition,  the  work  softening  feature  of  the  expansive  hardening 
function  may  he  modified  or  deleted  in  future  versions  of  the  conic  model 
to  insure  uniqueness  and  stability,  and  to  achieve  a  finite,  constant 
shear  strength  at  large  shear  strain  (a  non-zero  critical  state).  The 
model  is  called  a  conic  model  because  all  three  controlling  surfaces  in 
principal  stress  space  have  both  triaxial  and  octahedral  cross  sections 
which  are  conic  sections.  It  is  also  called  a  three  invariant  model 
because  the  expansive  yield  surface  involves  three  independent  stress 
invariants:  the  first  total  stress  invariant  and  the  second  and  third 
deviator  stress  invariants.  The  ARA  conic  model  rates  favorably  with 
respect  to  all  ten  evaluation  criteria  in  Table  3.1. 

At  present  the  conic  model  uses  the  incremental  stiffness  formulation 
developed  in  Appendix  G,  rather  than  a  trial  and  error  yield  surface 
violation  correction  procedure  such  as  that  discussed  in  Appendix  S  for 
the  cap  model.  However,  the  initial  strain  increments  needed  for 
numerical  stability  of  the  conic  model  are  very  small  (of  the  order  of 
10-^),  and  a  trial  and  error  correction  procedure  may  be  an  economic 
necessity  when  using  it  to  solve  dynamic  or  static  boundary  value  problems. 

W .?  Assumptions 

With  one  important  exception,  the  assumptions  underlying  the  conic 
model  are  the  same  as  those  underlying  the  Lade  model,  discussed  in 
Section  V.8.?.  The  one  exception  is  the  conic  model  enforces  the 
dissipation  condition,  whereas  the  Lade  model  does  not.  Using  two 
separate  yield  surfaces,  each  with  its  own  potential  surface  and  hardening 
rule,  assumes  a  material  really  has  two  separate  yield  surfaces.  The 
strain  hardening  parameter  for  each  yield  surface  is  the  corresponding 


plastic  work.  Of  course  all  these  elements  of  an  elastoplastic  model 
really  amount  to  physically  motivated  curve  fitting,  the  acid  tests  of 
which  are  predictive  accuracy  and  ease  of  use. 

W.:  Basic  Equations 

The  conic  model  compressive  yield  criterion  has  the  same  general  form 
as  does  the  Lade  model. 


fc  "  fc 


f"  =  0 
c 


but  the  stress  related  function,  f^,  has  a  slightly  different  form. 


i  2,-2  2 

'c  =  3o0CT  t0CT 


where 


(V.8.1) 

(W.l) 


r  =  compressive  yield  surface  ellipse  axis  ratio  (see  Figure  (U.l)). 
The  compressive  hardening  function,  f",  has  the  same  form 
as  does  the  Lade  model. 


where 


(V.8.3) 


Pa  =  atmospheric  pressure 
C,  p  =  material  parameters 

and  the  compressive  plastic  work,  Wc  is  defined  by  the  equation 

Wc  -  /HTM  fw'2) 

The  compressive  plastic  potential,  to  which  the  compressive  plastic  strain 
increment  vector,  |d e c],  is  normal,  is 

9c  =  fc  =  3o0CT2  +  3r2j  OCT  (W-3) 

The  conic  model  expansive  yield  criterion  has  the  same  general  form 

as  does  the  Lade  model. 


but  the  stress  related  function 


has  a  different  form. 


(V.8.6) 
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f  ■  =  (1  -  Ecos  3u)  —  +  m 

P  Vpa J  \  OCT  J 


(W.4) 


in  which 


f  1  _  f  I 

p  =  p,max  =  11 1 


at  failure 


(V.8.8) 


The  expansive  hardening  function,  f"  has  the  same  form  as  does  the  Lade 
mode  1 .  p 


-bW  /WnV/q 

irj 


( V.8.9) 


In  Equations  (W.4),  (V.8.8)  and  (V.8.9) 

E,  m,  n^.  a,  b,  q  =  material  parameters 


and  the  expansive  plastic  work,  W^,  is  defined  by  the  equation 


(  V.8.10) 


The  expansive  plastic  potential,  to  which  the  expansive  plastic  strain 
increment  vector,  {depJ  ,  is  normal,  is 


P  V  P, 


9„  =  |  -  Ecos  3«)  -  "V  ra 


* 


1  *  nC-M 


(W.5) 


where  is  also  a  material  parameter. 

The  unload! ng/rel oadi ng  elastic  modulus  is  assumed  to  be  given  by  the 
expression 

n 


Eur  =  Kurpai 


/°0CT\ 

\  pa  J 


(W.6! 


where 


Kur,  n  =  material  parameters 


A  V 


and  the  unloading/reloading  Poisson's  ratio, 


V 


jr  * 


is  computed  as 


described  in  Appendix  T,  or  assumed.  When  neither  yield  surface  is  active 
the  material  response  is  i ncremental ly  elastic. 

An  unloading/reloading  hysteresis  option  for  cyclic  loading  has  been 
formulated  but  not  yet  fully  implemented.  The  deviatoric  strain 
increments  are  considered  completely  elastic,  so  that  the  deviator  stress 
increments  are  given  by  the  equation 


ur 


ur 


de-  . 
U 


(W.7) 


However,  a  portion  of  the  volumetric  strain  increment  is  attributed  to 
unloadi ng/rel oadi ng  shear  hysteresis  rather  than  to  a  change  in  octahedral 
normal  stress.  In  other  words,  a  volumetric  strain  increment,  de,  would 
accompany  an  octahedral  shear  strain  increment,  dz,  even  under  constant 
octahedral  normal  stress.  If  the  octahedral  shear  strain  is 


z 


e  -  -e  •  . 

ij  U 


(W.8) 


then  the  octahedral  shear  strain  increment,  dz,  is  given  by  the  expressions 


dz 


e  .  -de  .  . 
U  iJ 

3z 


(z  >  0) 


(2  =  0) 


(W.9) 


The  octahedral  shear  strain,  z,  is  nonnegative,  but  the  octahedral  shear 
strain  increment,  dz,  can  be  positive,  zero,  or  negative.  The  volumetric 
strain  increment  due  to  shear,  de,  is  assumed  to  be  related  to  the 
octahedral  shear  strain  increment,  dz,  by  the  equation 


de  = 


dz 


(W.10) 


where 


A,  M,  x,  y,  e  =  material  parameters 

N  =  number  of  octahedral  shear  strain  reversals 
The  parameters  are  chosen  so  the  denominator  in  Equation  (W.10)  will 
always  be  positive.  Equation  (W.10)  can  accommodate  irregular  cyclic 
loading,  and  produces  hysteresis  loops  which  narrow  progressively  as  the 
number  of  octahedral  shear  strain  reversals,  N,  increases.  The  octahedral 
normal  stress  increment  is  calculated  by  the  equation 

do0CT  =  K(dE|(k  -  de)  (W.U) 

The  method  for  calculating  the  octahedral  cross  section  of  the  conic 
failure  surface  (the  expansive  yield  surface  at  its  maximum  extent)  is 
given  in  Appendix  U. 

The  computational  features  of  the  conic  model  are  covered  in  the 
appendices.  Using  the  value  of  Young's  modulus  from  Equation  (W.6)  and  an 
assumed  Poisson's  ratio,  the  elastic  incremental  stiffness  matrix  is  given 
by  Equation  (J.32).  The  general  equations  of  elastoplasticity  given  in 
Appendix  D  apply.  The  elastoplastic  incremental  stiffness  matrices  are 
calculated  as  described  in  Appendix  G,  and  the  polar  mode  check  described 
in  Appendix  I  is  used  when  the  stress  point  is  at  the  intersection  of  the 
two  yield  surfaces. 

W.fl  Parameter  Determination 

Determination  of  the  material  parameters  in  Equations  (V.8.1)  through 
(W.10)  is  discussed  in  Appendix  U.  Figures  (W.l)  through  (W.20) 
illustrate  the  fitting  process  for  the  ARA  model  using  remolded  CARES-DRY 
Sand  data.  All  the  fitting  was  performed  automatically  by  the  Soil 
Element  Model.  Figure  (W.l)  shows  the  variation  of  Young's  modulus  with 


maximum  past  pressure  for  this  material.  This  plot  represents  the  best 
straight  line  fit  to  observed  unload-reload  stiffness  from  isotropic 
compression,  uniaxial  strain,  and  triaxial  compression  tests  reported  by 
[Cargile  (1984)].  The  corresponding  pressures  are  the  test  pressures  at 
which  unloading  commenced. 

Figure  (W.?)  is  the  hydrostatic  data  (loading  only)  used  to  obtain 
the  compressive  plastic  work  parameters.  In  Figure  (W.3),  this  curve  has 
been  integrated  to  show  compressive  plastic  work  versus  pressure.  The 
integration  process  assumes  the  variation  of  elastic  properties  with 
confining  pressure  as  previously  determined.  The  data,  shown  as  open 
circles,  does  not  fall  along  a  single  straight  line,  even  in  log-log 
space.  Therefore,  several  straight  line  segments  have  been  fit,  as  shown, 
with  breakpoints  between  them  defined  on  the  basis  of  compressive  plastic 
work.  The  ARA  conic  model  is  presently  formulated  so  that  up  to  four 
segments  can  be  used. 

Figure  (W.4)  shows  triaxial  compression  stress-strain  data,  which,  at 
some  confining  pressures,  has  been  artificially  extended  to  near  peak 
response.  This  has  been  done  by  assuming  a  hyperbolic  shape  for  the  last 
twenty-five  percent  of  both  the  axial  and  radial  stress  difference-strain 
curves.  A“Kondner  plot  of  strain/stress  difference  vs.  strain  is  then 
used  to  extend  the  data  to  ninety-five  percent  of  peak  stress  difference 
as  predicted  by  the  straight  line  fit.  Two  such  plots  are  shown  in 
Figures  (W.5)  and  (W.6),  for  confining  pressures  of  3.4  MPa  and  7.0  MPa, 
respectively.  Some  type  of  data  extension  to  "failure"  is  necessary  for 
this  type  of  model  to  yield  reasonable  values  of  peak  expansive  work. 


w  m  1  ip  v 


Figure  (W.7)  defines  the  shape  and  maximum  extension  of  the  expansive 

I 

yield  surface,  f  Peak  values  have  been  plotted  for  triaxial  tests 
run  at  five  different  confining  pressures. 

Given  the  previously  determined  elastic  and  compressive  plastic  work 
behaviors.  Figures  (W.8)  and  (W.9)  show  the  predicted  compressive  and 
expansive  plastic  strains  during  the  triaxial  test.  Figures  (W.10) 
through  (W.14)  break  down  the  total  volumetric  strain  response  at  each 
confining  pressure  into  components.  Shown  are  the  total,  elastic,  and 
elastic  plus  compressive  plastic  strains.  Expansive  plastic  strain  is  the 
difference  between  the  total  and  elastic  plus  compressive  plastic  curves. 
According  to  the  model  theory,  this  should  always  be  a  negative  quantity, 
i.e.,  the  expansive  yield  surface  should  always  produce  expansive  volume 
strain  components.  At  some  confining  pressures,  however,  this  does  not 
hold  true.  [See  Figures  (W.ll),  (H.13)  and  This  implies  that 

the  volumetric  response  determined  from  the  isotropic  compression  test 
cannot  be  totally  representati ve  of  the  triaxial  test  volumetric 
response.  Note  that  this  can  be  remedied  by  using  an  envelope  of 
pressure-volume  behavior  which  encompasses  both  tests.  Further  adjustment 
of  the  r  factor  may  also  help. 

The  observed  shape  of  the  expansive  surface  hardening  function  is 
shown  in  Figure  (W.15)  for  each  confining  pressure.  Peak  values  from 
these  curves  are  plotted  in  Figure  (W.16),  and  shape-related  values  are 
plotted  in  Figure  (W.17). 

Parameters  for  the  expansive  plastic  potential  surface  are  derived 

ll 

from  a  plot  of  vs.  fp  for  several  confining  pressures,  as  shown 
in  Figure  (W.18).  The  data  are  derived  from  triaxial  volume  response  and 
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model  predictions,  and  a  plane  is  passed  through  the  entire  set.  The  very 
rough  appearance  of  the  data  is  due  to  three  factors: 


(a)  rough  hand  digitization  of  test  data  which  had  been  plotted  on  a 
small  area; 

(h)  elimination  of  unload-reload  loops  from  the  data,  which  did  in 
fact  have  a  significant  effect  on  volumetric  response  (the  large 
dips,  most  noticeable  at  g,  coincide  with  the  eliminated 
loops) ; 

(c)  previous  fitting  processes  which  may  introduce  some 

inconsistencies  between  test  data  and  predicted  model  response. 

Figures  (W.19)  and  (W.?0)  summarize  the  shape  of  the  fitted  expansive 
and  compressive  surfaces  in  the  triaxial  and  octahedral  planes, 
respectively.  Table  (W.l)  summarizes  the  ARA  parameters  for  remolded 


CARES-DRY  Sand. 

W.5  Calculated  Behavior 

The  multi -linear  fit  to  compressive  plastic  work  produces  an 
excellent  fit  to  loading  isotropic  compression  stress-strain  data,  as 
shown  in  Figure  (W.21).  Note  that  both  concave  and  convex  behavior  can  be 
matched  for  a  large  range  of  pressure.  The  model  currently  predicts  only 
linear  elastic  behavior  on  unloading-reloading,  so  the  observed  "tail"  on 
the  unload  curve  is  not  well  modeled. 

CTC  stress-strain  behavior  is  matched  fairly  accurately,  as  shown  in 
Figures  (W . 22 )  through  (W.25).  The  model  cannot  follow  the  uneven  volume 
strain  response  shown  in  Figure  (W.?A),  but  does  match  the  overall  trends 
of  compaction  and  subsequent  dilation  at  each  confining  pressure.  Figure 
•  (W.25)  shows  very  good  prediction  of  pressure- volume  behavior  for  this 

test. 
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Peak  values  of  stress  difference  for  the  CTE  test  were  used  for 
defining  the  asymmetry  of  the  expansive  yield  surface,  so  they  are  well 
matched,  as  shown  in  Figure  (W.25).  The  shape  of  the  stress-strain 
behavior  is  also  predicted  quite  well  [Figures  (W.26)  and  (W.27)]. 

Observed  volume  behavior  for  this  test  is  essentially  all  compaction,  and 
this  is  true  for  the  calculated  behavior  as  well  [Figures  (W.28)  and 
(W.29)].  Note,  however,  that  the  model  predicts  dilation  past  4-5  percent 
axial  strain,  while  the  samples  apparently  did  not  dilate  at  all. 

Calculated  RTC/RTE  behavior  [Figures  (W.30)  through  (W.33)] 
illustrates  several  important  features  of  the  ARA  model.  Stress-strain 
curves  [Figures  (W.30),  (W.31)  and  (W.33)]  are  smooth  for  both  tests. 
Substantial  softening  of  the  expansive  yield  surface  occurs  during  the  RTE 
tests.  As  shown  in  Figure  (W.31),  this  results  in  a  negative  shear 
stiffness  past  20-25  percent  strain  difference.  Dilation  is  predicted  for 
both  tests  [Figures  (W.32)  and  (W.33)]. 

PSC/PSE  calculations  show  non-symmetric  stress  strain  behavior 
[Figures  (W.3A)  and  (W.35)],  as  would  be  expected  due  to  the  non-symmetric 
yield  surface.  Note  that  no  work  softening  is  predicted,  because  W^  ^ 
has  not  been  achieved  for  any  of  the  confining  pressures.  Volume  strains 
are  predicted  to  be  initially  compressive  due  to  outward  movement  of  the 
elliptical  cap,  and  then  expansive  due  to  shear  dilation  [Figure  (W.36)]. 
Figure  (W.37)  summarizes  calculated  behavior  of  the  ARA  model  for  several 
different  tests  run  in  the  triaxial  device.  All  tests  start  at  7  MPa 
initial  confining  pressure. 

Uniaxial  strain  (UXC)  data  vs.  calculation  comparisons  are  shown  in 
Figures  (W.38)  through  (W.40).  Stress-strain  behavior  is  only  marginally 


good,  because  of  the  model's  tendency  to  stiffen  (due  to  dilation)  under 
non-isotropic  stress  paths.  The  stress  path  is  predicted  quite  well,  with 
the  exception  of  the  non-linear  unload-reload  excursions. 

The  observed  shape  of  the  UXE  stress  path  is  predicted  very  well  by 
the  ARA  model  as  shown  in  Figure  (W.41).  Apparently,  the  ultimate 
expansive  yiel d  surface  position  is  somewhat  more  expanded  than  predicted 
by  the  model.  Figure  (W.42)  shows  calculated  UXE  stress-strain  response. 

Axisymmetric  strain  path  results  for  the  ARA  model  [Figures  (W.43) 
through  (W.^F)]  indicate  that  the  determined  parameters  produce  behavior 
which  is  too  stiff  at  the  WES  level  of  confining  pressure  (7  MPa).  This 
point  can  be  confirmed  by  re-comparing  UXC  calculation  and  test  data  at  an 
axial  stress  of  about  12  MPa  (corresponding  to  a  pressure  of  about  7  MPa) 
in  Figure  (W.3P).  The  stresses  from  the  WES  strain  paths  [Figure  (W.43)] 
are  substantially  overpredicted  as  a  result.  Stress  difference  also  drops 
too  quickly,  as  radial  expansion  is  intensified.  The  model  does  somewhat 
better  with  the  Lade  strain  paths  [Figures  (W.44)  and  (W.45)],  perhaps 
because  the  initial  confining  pressure  is  lower.  In  fact,  for  these 
paths,  both  shear  and  volume  stiffness  are  somewhat  underpredicted. 

Figure  (W.*6)  shows  the  predicted  stress  path  for  the  true  triaxial 
test  strain  path,  and  Figure  (W.47)  shows  the  predicted  pressure- vol ume 


response. 


TABLE  W.l(a)  ARA  MODEL  PARAMETERS  FOR 
REMOLDED  CARES-DRY  SAND 


Parameter 

Symbol 

Variable 

Value 

Units 

Modulus  Constant 

^ur 

AKUR 

36.35 

Modulus  Exponent 

n 

AH 

0.8412 

-- 

Poisson's  Ratio 

ve 

AP01 

0.200 

— 

No.  fc  Segments 

ncrv 

ACRV 

3 

-- 

Hardening  Constant 

Cl 

AACC(l) 

4.645x10-5 

Hardening  Exponent  Seg.  1 

Pi 

AAPC(l) 

1.401 

-- 

Hardening  Breakpoint 

hk2 

ABRK{1) 

2.195 

-- 

c2 

AACC ( 2 ) 

6.086xl0-? 

Seg.  ? 

P2 

AAPCI2) 

0.4667 

-- 

bk2 

AB  RK ( 2 ) 

22.84 

-- 

C3 

AACC  ( 3 ) 

7.516x10-1 

-- 

Seg.  3 

P3 

AA PC(3) 

0.2688 

— 

CAP  Axes  Ratio 

r 

AR 

0.250 

-- 

Yield  Constant 

E 

AEY 

0.111 

-- 

Yield  Exponent 

m 

AMY 

2.875x10-4 

-- 

Failure  Constant 

n 

AETA1 

0.6454 

-- 

P 

APBAR 

0.5057 

_  _ 

Work  Hardening  Constants 

1 

AL 

0.8691 

-- 

a 

AALPH 

5.000 

-- 

8 

ABETA 

-2.631x10-3 

-- 

t 

ATG 

-0.9646 

— 

Plastic  Potential  Constants 

R 

ARG 

2.182x10-3 

S 

ASG 

1.860 

-- 

Plastic  Expansive  Work 

°0CT,max 

ASOCT 

Pa 

at  *p,max 

Wp,pk 

AWPPK 

Pa 

Work  Hardening  Constants 

q 

AO 

varies,  see 

-- 

a 

AA 

Table  W.l(b) 

-- 

b 

AB 

-- 

Plastic  Potential  Constant 

^2,0 

AETA20 

-  _ 

Initial  Plastic  Comp.  Work 

WC,0 

AWC 

Pa 

Mass  Density 

p 

RHOREF 

1900 

kg/m3 

TABLE  W.l(b)  CONFINING  PRESSURE  DEPENDENT 
ARA  MODEL  PARAMETERS 


°3c 

u 

p,pk 

q 

a 

b 

n2,0 

Wc,0 

(MPa)  (Pa) 

(Pa) 

A 

P 

0.1 

5.07xl04 

4.997 

0.9058 

3.950x10  * 

-0.962a 

3.03x10 

0.4 

1.74xl05 

4.989 

0.7076 

1.151xl0"6 

-0.9602 

1.62x10 

1.8 

6.25xl05 

4.953 

0.5471 

3. 232xl0'7 

-0.9554 

9.97x10 

3.5 

l.llxl O6 

4.909 

0.4856 

1.830xl0"7 

-0.9518 

6.42x10 

7.0 

2. 03xl06 

4.818 

0.4263 

1.021xl0"7 

-0.9465 

6.87x10 

20.0 

5. 06x1 06 

4.481 

0.3370 

4.407x10"® 

-0.9340 

1.21x10 

32.0 

7.62xl06 

4.170 

0.2911 

3.148x10'® 

-0.9258 

1.56x10 

59.0 

1.29xl07 

3.47a 

0.2132 

2.227x10'® 

-0.9121 

2.16x10 

100.0 

2. OFxlO7 

2.405 

0.1075 

2.027x10'® 

-0.8961 

2.87x10 
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